Mapping transiently formed and sparsely populated conformations on a complex energy landscape

Impact of surface-active ionic solutions on the structure and function of laccase from trametes versicolor: Insights from molecular dynamics simulations

Many protein-ionic liquid investigations have examined laccase interactions. Laccases are a class of poly-copper oxidoreductases that retain significant biotechnological relevance owing to their notable oxidative capabilities and their application in the elimination of synthetic dyes, phenolic compounds, insecticides, and various other substances. This study investigates the impact of surface active ionic liquids (SAILs), namely, decyltrimethylammonium bromide [N10111][Br] and 1-decyl-3-methylimidazolium chloride [C10mim][Cl] as cationic surfactant ionic liquids and cholinium decanoate [Chl][Dec], an anionic surfactant ionic liquid, on the structure and function of laccase from the fungus Trametes versicolor (TvL) by the molecular dynamics (MD) simulation method. In summary, this study showed that laccase solvent-accessible surface area increased in the ionic liquid [Chl][Dec] while it decreased in the other two ionic liquids. Interestingly, [Chl][Dec] ionic liquid components formed hydrogen bonds with laccase, while [N10111][Br] and [C10mim][Cl] components were unable to form hydrogen bonds with laccase. The quantity of hydrogen bonds formed between water molecules and the enzyme was also diminished in the presence of [Chl][Dec] in comparison to the other two ionic liquids. especially at a concentration of 250 mM. In 250 mM concentrations of [N10111][Br] and [C10mim][Cl], clusters of long-chain cations are likely to form near the copper T1 site. However, even at low [Chl][Dec] concentrations, long [Dec]- chains were observed to penetrate the enzyme near the copper T1 site, and at 250 mM [Chl][Dec], a large cluster of anions occupied the opening of the active site. The results of the analysis also show that the interaction between the [Dec]- anion and the enzyme is stronger than the interaction between [N10111]+ and [C10mim]+ with laccase; in addition, the [Dec]- anion, compared to [Br]- and [Cl]- has a much greater tendency to bind with the enzyme residues.

Structural dynamics in chromatin unraveling by pioneer transcription factors

Pioneer transcription factors are proteins with a dual function. First, they regulate transcription by binding to nucleosome-free DNA regulatory elements. Second, they bind to DNA while wrapped around histone proteins in the chromatin and mediate chromatin opening. The molecular mechanisms that connect the two functions are yet to be discovered. In recent years, pioneer factors received increased attention mainly because of their crucial role in promoting cell fate transitions that could be used for regenerative therapies. For example, the three factors required to induce pluripotency in somatic cells, Oct4, Sox2, and Klf4 were classified as pioneer factors and studied extensively. With this increased attention, several structures of complexes between pioneer factors and chromatin structural units (nucleosomes) have been resolved experimentally. Furthermore, experimental and computational approaches have been designed to study two unresolved, key scientific questions: First, do pioneer factors induce directly local opening of nucleosomes and chromatin fibers upon binding? And second, how do the unstructured tails of the histones impact the structural dynamics involved in such conformational transitions? Here we review the current knowledge about transcription factor-induced nucleosome dynamics and the role of the histone tails in this process. We discuss what is needed to bridge the gap between the static views obtained from the experimental structures and the key structural dynamic events in chromatin opening. Finally, we propose that integrating nuclear magnetic resonance spectroscopy with molecular dynamics simulations is a powerful approach to studying pioneer factor-mediated dynamics of nucleosomes and perhaps small chromatin fibers using native DNA sequences.

Structural basis for sugar perception by Drosophila gustatory receptors

Insects rely on a family of seven transmembrane proteins called gustatory receptors (GRs) to encode different taste modalities, such as sweet and bitter. We report structures of Drosophila sweet taste receptors GR43a and GR64a in the apo and sugar-bound states. Both GRs form tetrameric sugar-gated cation channels composed of one central pore domain (PD) and four peripheral ligand-binding domains (LBDs). Whereas GR43a is specifically activated by the monosaccharide fructose that binds to a narrow pocket in LBDs, disaccharides sucrose and maltose selectively activate GR64a by binding to a larger and flatter pocket in LBDs. Sugar binding to LBDs induces local conformational changes, which are subsequently transferred to the PD to cause channel opening. Our studies reveal a structural basis for sugar recognition and activation of GRs.

Kinetics from Metadynamics: Principles, Applications, and Outlook

Metadynamics is a popular enhanced sampling algorithm for computing the free energy landscape of rare events by using molecular dynamics simulation. Ten years ago, Tiwary and Parrinello introduced the infrequent metadynamics approach for calculating the kinetics of transitions across free energy barriers. Since then, metadynamics-based methods for obtaining rate constants have attracted significant attention in computational molecular science. Such methods have been applied to study a wide range of problems, including protein-ligand binding, protein folding, conformational transitions, chemical reactions, catalysis, and nucleation. Here, we review the principles of elucidating kinetics from metadynamics-like approaches, subsequent methodological developments in this area, and successful applications on chemical, biological, and material systems. We also highlight the challenges of reconstructing accurate kinetics from enhanced sampling simulations and the scope of future developments.

Uncovering the Kinetic Characteristics and Degradation Preference of PROTAC Systems with Advanced Theoretical Analyses

Proteolysis-targeting chimeras (PROTACs), which can selectively induce the degradation of target proteins, represent an attractive technology in drug discovery. A large number of PROTACs have been reported, but due to the complicated structural and kinetic characteristics of the target-PROTAC-E3 ligase ternary interaction process, the rational design of PROTACs is still quite challenging. Here, we characterized and analyzed the kinetic mechanism of MZ1, a PROTAC that targets the bromodomain (BD) of the bromodomain and extra terminal (BET) protein (Brd2, Brd3, or Brd4) and von Hippel-Lindau E3 ligase (VHL), from the kinetic and thermodynamic perspectives of view by using enhanced sampling simulations and free energy calculations. The simulations yielded satisfactory predictions on the relative residence time and standard binding free energy (r_p > 0.9) for MZ1 in different Brd^BD-MZ1-VHL ternary complexes. Interestingly, the simulation of the PROTAC ternary complex disintegration illustrates that MZ1 tends to remain on the surface of VHL with the BD proteins dissociating alone without a specific dissociation direction, indicating that the PROTAC prefers more to bind with E3 ligase at the first step in the formation of the target-PROTAC-E3 ligase ternary complex. Further exploration of the degradation difference of MZ1 in different Brd systems shows that the PROTAC with higher degradation efficiency tends to leave more lysine exposed on the target protein, which is guaranteed by the stability (binding affinity) and durability (residence time) of the target-PROTAC-E3 ligase ternary complex. It is quite possible that the underlying binding characteristics of the Brd^BD-MZ1-VHL systems revealed by this study may be shared by different PROTAC systems as a general rule, which may accelerate rational PROTAC design with higher degradation efficiency.

Systematic conformation-to-phenotype mapping via limited deep sequencing of proteins

Non-native conformations drive protein-misfolding diseases, complicate bioengineering efforts, and fuel molecular evolution. No current experimental technique is well suited for elucidating them and their phenotypic effects. Especially intractable are the transient conformations populated by intrinsically disordered proteins. We describe an approach to systematically discover, stabilize, and purify native and non-native conformations, generated in vitro or in vivo, and directly link conformations to molecular, organismal, or evolutionary phenotypes. This approach involves high-throughput disulfide scanning (HTDS) of the entire protein. To reveal which disulfides trap which chromatographically resolvable conformers, we devised a deep-sequencing method for double-Cys variant libraries of proteins that precisely and simultaneously locates both Cys residues within each polypeptide. HTDS of the abundant E. coli periplasmic chaperone HdeA revealed distinct classes of disordered hydrophobic conformers with variable cytotoxicity depending on where the backbone was cross-linked. HTDS can bridge conformational and phenotypic landscapes for many proteins that function in disulfide-permissive environments.

Insights into membrane association of the SMP domain of extended synaptotagmin

The Synaptotagmin-like Mitochondrial-lipid-binding Protein (SMP) domain is a newly identified lipid transfer module present in proteins that regulate lipid homeostasis at membrane contact sites (MCSs). However, how the SMP domain associates with the membrane to extract and unload lipids is unclear. Here, we performed in vitro DNA brick-assisted lipid transfer assays and in silico molecular dynamics simulations to investigate the molecular basis of the membrane association by the SMP domain of extended synaptotagmin (E-Syt), which tethers the tubular endoplasmic reticulum (ER) to the plasma membrane (PM). We demonstrate that the SMP domain uses its tip region to recognize the extremely curved subdomain of tubular ER and the acidic-lipid-enriched PM for highly efficient lipid transfer. Supporting these findings, disruption of these mechanisms results in a defect in autophagosome biogenesis contributed by E-Syt. Our results suggest a model that provides a coherent picture of the action of the SMP domain at MCSs.

Ligand Gaussian Accelerated Molecular Dynamics 2 (LiGaMD2): Improved Calculations of Ligand Binding Thermodynamics and Kinetics with Closed Protein Pocket

Ligand binding thermodynamics and kinetics are critical parameters for drug design. However, it has proven challenging to efficiently predict ligand binding thermodynamics and kinetics from molecular simulations due to limited simulation timescales. Protein dynamics, especially in the ligand binding pocket, often plays an important role in ligand binding. Based on our previously developed Ligand Gaussian accelerated molecular dynamics (LiGaMD), here we present LiGaMD2 in which a selective boost potential was applied to both the ligand and protein residues in the binding pocket to improve sampling of ligand binding and dissociation. To validate the performance of LiGaMD2, the T4 lysozyme (T4L) mutants with open and closed pockets bound by different ligands were chosen as model systems. LiGaMD2 could efficiently capture repetitive ligand dissociation and binding within microsecond simulations of all T4L systems. The obtained ligand binding kinetic rates and free energies agreed well with available experimental values and previous modeling results. Therefore, LiGaMD2 provides an improved approach to sample opening of closed protein pockets for ligand dissociation and binding, thereby allowing for efficient calculations of ligand binding thermodynamics and kinetics.

Systematic conformation-to-phenotype mapping via limited deep-sequencing of proteins

Non-native conformations drive protein misfolding diseases, complicate bioengineering efforts, and fuel molecular evolution. No current experimental technique is well-suited for elucidating them and their phenotypic effects. Especially intractable are the transient conformations populated by intrinsically disordered proteins. We describe an approach to systematically discover, stabilize, and purify native and non-native conformations, generated in vitro or in vivo, and directly link conformations to molecular, organismal, or evolutionary phenotypes. This approach involves high-throughput disulfide scanning (HTDS) of the entire protein. To reveal which disulfides trap which chromatographically resolvable conformers, we devised a deep-sequencing method for double-Cys variant libraries of proteins that precisely and simultaneously locates both Cys residues within each polypeptide. HTDS of the abundant E. coli periplasmic chaperone HdeA revealed distinct classes of disordered hydrophobic conformers with variable cytotoxicity depending on where the backbone was cross-linked. HTDS can bridge conformational and phenotypic landscapes for many proteins that function in disulfide-permissive environments.

Automated Path Searching Reveals the Mechanism of Hydrolysis Enhancement by T4 Lysozyme Mutants

Bacteriophage T4 lysozyme (T4L) is a glycosidase that is widely applied as a natural antimicrobial agent in the food industry. Due to its wide applications and small size, T4L has been regarded as a model system for understanding protein dynamics and for large-scale protein engineering. Through structural insights from the single conformation of T4L, a series of mutations (L99A,G113A,R119P) have been introduced, which have successfully raised the fractional population of its only hydrolysis-competent excited state to 96%. However, the actual impact of these substitutions on its dynamics remains unclear, largely due to the lack of highly efficient sampling algorithms. Here, using our recently developed travelling-salesman-based automated path searching (TAPS), we located the minimum-free-energy path (MFEP) for the transition of three T4L mutants from their ground states to their excited states. All three mutants share a three-step transition: the flipping of F114, the rearrangement of α0/α1 helices, and final refinement. Remarkably, the MFEP revealed that the effects of the mutations are drastically beyond the expectations of their original design: (a) the G113A substitution not only enhances helicity but also fills the hydrophobic Cavity I and reduces the free energy barrier for flipping F114; (b) R119P barely changes the stability of the ground state but stabilizes the excited state through rarely reported polar contacts S117_OG:N132_ND2, E11_OE1:R145_NH1, and E11_OE2:Q105_NE2; (c) the residue W138 flips into Cavity I and further stabilizes the excited state for the triple mutant L99A,G113A,R119P. These novel insights that were unexpected in the original mutant design indicated the necessity of incorporating path searching into the workflow of rational protein engineering.

Transient exposure of a buried phosphorylation site in an autoinhibited protein

Autoinhibition is a mechanism used to regulate protein function, often by making functional sites inaccessible through the interaction with a cis-acting inhibitory domain. Such autoinhibitory domains often display a substantial degree of structural disorder when unbound, and only become structured in the inhibited state. These conformational dynamics make it difficult to study the structural origin of regulation, including effects of regulatory post-translational modifications. Here, we study the autoinhibition of the Dbl Homology domain in the protein Vav1 by the so-called acidic inhibitory domain. We use molecular simulations to study the process by which a mostly unstructured inhibitory domain folds upon binding and how transient exposure of a key buried tyrosine residue makes it accessible for phosphorylation. We show that the inhibitory domain, which forms a helix in the bound and inhibited stated, samples helical structures already before binding and that binding occurs via a molten-globule-like intermediate state. Together, our results shed light on key interactions that enable the inhibitory domain to sample a finely tuned equilibrium between an inhibited and a kinase-accessible state.

Energy penalties enhance flexible receptor docking in a model cavity

Protein flexibility remains a major challenge in library docking because of difficulties in sampling conformational ensembles with accurate probabilities. Here, we use the model cavity site of T4 lysozyme L99A to test flexible receptor docking with energy penalties from molecular dynamics (MD) simulations. Crystallography with larger and smaller ligands indicates that this cavity can adopt three major conformations: open, intermediate, and closed. Since smaller ligands typically bind better to the cavity site, we anticipate an energy penalty for the cavity opening. To estimate its magnitude, we calculate conformational preferences from MD simulations. We find that including a penalty term is essential for retrospective ligand enrichment; otherwise, high-energy states dominate the docking. We then prospectively docked a library of over 900,000 compounds for new molecules binding to each conformational state. Absent a penalty term, the open conformation dominated the docking results; inclusion of this term led to a balanced sampling of ligands against each state. High ranked molecules were experimentally tested by Tm upshift and X-ray crystallography. From 33 selected molecules, we identified 18 ligands and determined 13 crystal structures. Most interesting were those bound to the open cavity, where the buried site opens to bulk solvent. Here, highly unusual ligands for this cavity had been predicted, including large ligands with polar tails; these were confirmed both by binding and by crystallography. In docking, incorporating protein flexibility with thermodynamic weightings may thus access new ligand chemotypes. The MD approach to accessing and, crucially, weighting such alternative states may find general applicability.

Ligand unbinding mechanisms and kinetics for T4 lysozyme mutants from τRAMD simulations

The protein-ligand residence time, τ, influences molecular function in biological networks and has been recognized as an important determinant of drug efficacy. To predict τ, computational methods must overcome the problem that τ often exceeds the timescales accessible to conventional molecular dynamics (MD) simulation. Here, we apply the τ-Random Acceleration Molecular Dynamics (τRAMD) method to a set of kinetically characterized complexes of T4 lysozyme mutants with small, engineered binding cavities. τRAMD yields relative ligand dissociation rates in good accordance with experiments across this diverse set of complexes that differ with regard to measurement temperature, ligand identity, protein mutation and binding cavity. τRAMD thereby allows a comprehensive characterization of the ligand egress routes and determinants of τ. Although ligand dissociation by multiple egress routes is observed, we find that egress via the predominant route determines the value of τ. We also find that the presence of a greater number of metastable states along egress pathways leads to slower protein-ligand dissociation. These physical insights could be exploited in the rational optimization of the kinetic properties of drug candidates.

•

Relative residence times are computed for T4 lysozyme mutant-ligand complexes.

•

τ-Random Acceleration Molecular Dynamics provide efficient sampling of unbinding.

•

Computed dissociation rates show good agreement with all available measured values.

•

Ligand egress via the predominant route determines the value of the residence time.

•

The presence of metastable states along egress pathways slows down dissociation.

DEER-PREdict: Software for efficient calculation of spin-labeling EPR and NMR data from conformational ensembles

Owing to their plasticity, intrinsically disordered and multidomain proteins require descriptions based on multiple conformations, thus calling for techniques and analysis tools that are capable of dealing with conformational ensembles rather than a single protein structure. Here, we introduce DEER-PREdict, a software program to predict Double Electron-Electron Resonance distance distributions as well as Paramagnetic Relaxation Enhancement rates from ensembles of protein conformations. DEER-PREdict uses an established rotamer library approach to describe the paramagnetic probes which are bound covalently to the protein.DEER-PREdict has been designed to operate efficiently on large conformational ensembles, such as those generated by molecular dynamics simulation, to facilitate the validation or refinement of molecular models as well as the interpretation of experimental data. The performance and accuracy of the software is demonstrated with experimentally characterized protein systems: HIV-1 protease, T4 Lysozyme and Acyl-CoA-binding protein. DEER-PREdict is open source (GPLv3) and available at github.com/KULL-Centre/DEERpredict and as a Python PyPI package pypi.org/project/DEERPREdict.

Investigating Conformational Dynamics and Allostery in the p53 DNA-Binding Domain Using Molecular Simulations

The p53 tumor suppressor is a multifaceted context-dependent protein, which is involved in multiple cellular pathways, with the ability to either keep the cells alive or to kill them through mechanisms such as apoptosis. To complicate this picture, cancer cells that express mutant p53 becomes addicted to the mutant activity, so that the mutant variant features a myriad of gain-of-function activities, opening different venues for therapy. This makes essential to think outside the box and apply new approaches to the study of p53 structure-(mis)function relationship to find new critical components of its pathway or to understand how known parts are interconnected, compete, or cooperate. In this context, I will here illustrate how to integrate different computational methods to the identification of possible allosteric effects transmitted from the DNA binding interface of p53 to regions for cofactor recruitment. The protocol can be extended to any other cases of study. Indeed, it does not necessarily apply only to the study of DNA-induced effects, but more broadly to the investigation of long-range effects induced by a biological partner that binds to a biomolecule of interest.

Combining Virtual Reality Visualization with Ensemble Molecular Dynamics to Study Complex Protein Conformational Changes

Molecular dynamics (MD) simulations are increasingly used to elucidate relationships between protein structure, dynamics, and their biological function. Currently, it is extremely challenging to perform MD simulations of large-scale structural rearrangements in proteins that occur on millisecond timescales or beyond, as this requires very significant computational resources, or the use of cumbersome "collective variable" enhanced sampling protocols. Here, we describe a framework that combines ensemble MD simulations and virtual reality visualization (eMD-VR) to enable users to interactively generate realistic descriptions of large amplitude, millisecond timescale protein conformational changes in proteins. Detailed tests demonstrate that eMD-VR substantially decreases the computational cost of folding simulations of a WW domain, without the need to define collective variables a priori. We further show that eMD-VR generated pathways can be combined with Markov state models to describe the thermodynamics and kinetics of large-scale loop motions in the enzyme cyclophilin A. Our results suggest eMD-VR is a powerful tool for exploring protein energy landscapes in bioengineering efforts.

Wepy: A Flexible Software Framework for Simulating Rare Events with Weighted Ensemble Resampling

Here, we introduce the open-source software framework wepy (https://github.com/ADicksonLab/wepy) which is a toolkit for running and analyzing weighted ensemble (WE) simulations. The wepy toolkit is in pure Python and as such is highly portable and extensible, making it an excellent platform to develop and use new WE resampling algorithms such as WExplore, REVO, and others while leveraging the entire Python ecosystem. In addition, wepy simplifies WE-specific analyses by defining out-of-core tree-like data structures using the cross-platform HDF5 file format. In this paper, we discuss the motivations and challenges for simulating rare events in biomolecular systems. As has previously been shown, high-dimensional WE resampling algorithms such as WExplore and REVO have been successful at these tasks, especially for rare events that are difficult to describe by one or two collective variables. We explain in detail how wepy facilitates implementation of these algorithms, as well as aids in analyzing the unique structure of WE simulation results. To explain how wepy and WE work in general, we describe the mathematical formalism of WE, an overview of the architecture of wepy, and provide code examples of how to construct, run, and analyze simulation results for a protein-ligand system (T4 Lysozyme in an implicit solvent). This paper is written with a variety of readers in mind, including (1) those curious about how to leverage WE rare-event simulations for their domain, (2) current WE users who want to begin using new high-dimensional resamplers such as WExplore and REVO, and (3) expert users who would like to prototype or implement their own algorithms that can be easily adopted by others.

Protein Dynamics Enables Phosphorylation of Buried Residues in Cdk2/Cyclin-A-Bound p27

Proteins carry out a wide range of functions that are tightly regulated in space and time. Protein phosphorylation is the most common post-translation modification of proteins and plays a key role in the regulation of many biological processes. The finding that many phosphorylated residues are not solvent exposed in the unphosphorylated state opens several questions for understanding the mechanism that underlies phosphorylation and how phosphorylation may affect protein structures. First, because kinases need access to the phosphorylated residue, how do such buried residues become modified? Second, once phosphorylated, what are the structural effects of phosphorylation of buried residues, and do they lead to changed conformational dynamics? We have used the ternary complex between p27^Kip1 (p27), Cdk2, and cyclin A to study these questions using enhanced sampling molecular dynamics simulations. In line with previous NMR and single-molecule fluorescence experiments, we observe transient exposure of Tyr88 in p27, even in its unphosphorylated state. Once Tyr88 is phosphorylated, we observe a coupling to a second site, thus making Tyr74 more easily exposed and thereby the target for a second phosphorylation step. Our observations provide atomic details on how protein dynamics plays a role in modulating multisite phosphorylation in p27, thus supplementing previous experimental observations. More generally, we discuss how the observed phenomenon of transient exposure of buried residues may play a more general role in regulating protein function.

Targeting a critical step in fungal hexosamine biosynthesis

Aspergillus fumigatus is a human opportunistic fungal pathogen whose cell wall protects it from the extracellular environment including host defenses. Chitin, an essential component of the fungal cell wall, is synthesized from UDP-GlcNAc produced in the hexosamine biosynthetic pathway. As this pathway is critical for fungal cell wall integrity, the hexosamine biosynthesis enzymes represent potential targets of antifungal drugs. Here, we provide genetic and chemical evidence that glucosamine 6-phosphate N-acetyltransferase (Gna1), a key enzyme in this pathway, is an exploitable antifungal drug target. GNA1 deletion resulted in loss of fungal viability and disruption of the cell wall, phenotypes that could be rescued by exogenous GlcNAc, the product of the Gna1 enzyme. In a murine model of aspergillosis, the Δgna1 mutant strain exhibited attenuated virulence. Using a fragment-based approach, we discovered a small heterocyclic scaffold that binds proximal to the Gna1 active site and can be optimized to a selective submicromolar binder. Taken together, we have provided genetic, structural, and chemical evidence that Gna1 is an antifungal target in A. fumigatus.

Local frustration determines loop opening during the catalytic cycle of an oxidoreductase

Local structural frustration, the existence of mutually exclusive competing interactions, may explain why some proteins are dynamic while others are rigid. Frustration is thought to underpin biomolecular recognition and the flexibility of protein-binding sites. Here, we show how a small chemical modification, the oxidation of two cysteine thiols to a disulfide bond, during the catalytic cycle of the N-terminal domain of the key bacterial oxidoreductase DsbD (nDsbD), introduces frustration ultimately influencing protein function. In oxidized nDsbD, local frustration disrupts the packing of the protective cap-loop region against the active site allowing loop opening. By contrast, in reduced nDsbD the cap loop is rigid, always protecting the active-site thiols from the oxidizing environment of the periplasm. Our results point toward an intricate coupling between the dynamics of the active-site cysteines and of the cap loop which modulates the association reactions of nDsbD with its partners resulting in optimized protein function.

Perturb-Scan-Pull: A Novel Method Facilitating Conformational Transitions in Proteins

Conformational transitions in proteins facilitate precise physiological functions. Therefore, it is crucial to understand the mechanisms underlying these processes to modulate protein function. Yet, studying structural and dynamical properties of proteins is notoriously challenging due to the complexity of the underlying potential energy surfaces (PES). We have previously developed the perturbation-response scanning (PRS) method to identify key residues that participate in the communication network responsible for specific conformational transitions. PRS is based on a residue-by-residue scan of the protein to determine the subset of residues/forces which provide the closest conformational change leading to a target conformational state, inasmuch as linear response theory applies to these motions. Here, we develop a novel method to further evaluate if conformational transitions may be triggered on the PES. We aim to study functionally relevant conformational transitions in proteins by using results obtained from PRS and feeding them as inputs to steered molecular dynamics simulations. The success and the transferability of the method are evaluated on three protein systems having different complexities of motion on the PES: calmodulin, adenylate kinase, and bacterial ferric binding protein. We find that the method captures the target conformation, while providing key residues and the optimum paths with relatively low free energy profiles.

Dynamic design: manipulation of millisecond timescale motions on the energy landscape of cyclophilin A

Proteins need to interconvert between many conformations in order to function, many of which are formed transiently, and sparsely populated. Particularly when the lifetimes of these states approach the millisecond timescale, identifying the relevant structures and the mechanism by which they interconvert remains a tremendous challenge. Here we introduce a novel combination of accelerated MD (aMD) simulations and Markov state modelling (MSM) to explore these 'excited' conformational states. Applying this to the highly dynamic protein CypA, a protein involved in immune response and associated with HIV infection, we identify five principally populated conformational states and the atomistic mechanism by which they interconvert. A rational design strategy predicted that the mutant D66A should stabilise the minor conformations and substantially alter the dynamics, whereas the similar mutant H70A should leave the landscape broadly unchanged. These predictions are confirmed using CPMG and R1ρ solution state NMR measurements. By efficiently exploring functionally relevant, but sparsely populated conformations with millisecond lifetimes in silico, our aMD/MSM method has tremendous promise for the design of dynamic protein free energy landscapes for both protein engineering and drug discovery.

The investigation of the G-quadruplex aptamer selectivity to Pb2+ ion: a joint molecular dynamics simulation and density functional theory study

The aptamers with the ability to form a G-quadruplex structure can be stable in the presence of some ions. Hence, study of the interactions between such aptamers and ions can be beneficial to determine the highest selective aptamer toward an ion. In this article, molecular dynamics (MD) simulations and quantum mechanics (QM) calculations have been applied to investigate the selectivity of the T30695 aptamer toward Pb²⁺ in comparison with some ions. The Free Energy Landscape (FEL) analysis indicates that Pb²⁺ has remained inside the aptamer during the MD simulation, while the other ions have left it. The Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) binding energies prove that the conformational stability of the aptamer is the highest in the presence of Pb²⁺. According to the compaction parameters, the greatest compressed ion-aptamer complex, and hence, the highest ion-aptamer interaction have been induced in the presence of Pb²⁺. The contact maps clarify the closer contacts between the nucleotides of the aptamer in the presence of Pb²⁺. The density functional theory (DFT) results show that Pb²⁺ forms the most stable complex with the aptamer, which is consistent with the MD results. The QM calculations reveal that the N-H bonds and the O…H distances are the longest and the shortest, respectively, in the presence of Pb²⁺. The obtained results verify that the strongest hydrogen bonds (HBs), and hence, the most compressed aptamer structure are induced by Pb²⁺. Besides, atoms in molecules (AIM) and natural bond orbital (NBO) analyses confirm the results.Communicated by Ramaswamy H. Sarma.

Metadynamics to Enhance Sampling in Biomolecular Simulations

Molecular dynamics is a powerful simulation method to provide detailed atomic-scale insight into a range of biological processes including protein folding, biochemical reactions, ligand binding, and many others. Over the last several decades, enhanced sampling methods have been developed to address the large separation in time scales between a molecular dynamics simulation (usually microseconds or shorter) and the time scales of biological processes (often orders of magnitude longer). This chapter specifically focuses on the metadynamics family of methods, which achieves enhanced sampling through the introduction of a history-dependent bias potential that is based on one or more slow degrees of freedom, called collective variables. We introduce the method and its recent variants related to biomolecular studies and then discuss frontier areas of the method. A large part of this chapter is devoted to helping new users of the method understand how to choose metadynamics parameters properly and apply the method to their system of interest.

Escape of a Small Molecule from Inside T4 Lysozyme by Multiple Pathways

The T4 lysozyme L99A mutant is often used as a model system to study small-molecule binding to proteins, but pathways for ligand entry and exit from the buried binding site and the associated protein conformational changes have not been fully resolved. Here, molecular dynamics simulations were employed to model benzene exit from its binding cavity using the weighted ensemble (WE) approach to enhance sampling of low-probability unbinding trajectories. Independent WE simulations revealed four pathways for benzene exit, which correspond to transient tunnels spontaneously formed in previous simulations of apo T4 lysozyme. Thus, benzene unbinding occurs through multiple pathways partially created by intrinsic protein structural fluctuations. Motions of several α-helices and side chains were involved in ligand escape from metastable microstates. WE simulations also provided preliminary estimates of rate constants for each exit pathway. These results complement previous works and provide a semiquantitative characterization of pathway heterogeneity for binding of small molecules to proteins.

Mechanisms for Benzene Dissociation through the Excited State of T4 Lysozyme L99A Mutant

The atomic-level mechanisms that coordinate ligand release from protein pockets are only known for a handful of proteins. Here, we report results from accelerated molecular dynamics simulations for benzene dissociation from the buried cavity of the T4 lysozyme Leu99Ala mutant (L99A). In these simulations, benzene is released through a previously characterized, sparsely populated room-temperature excited state of the mutant, explaining the coincidence for experimentally measured benzene off rate and apo protein slow-timescale NMR relaxation rates between ground and excited states. The path observed for benzene egress is a multistep ligand migration from the buried cavity to ultimate release through an opening between the F/G-, H-, and I-helices and requires a number of cooperative multiresidue and secondary-structure rearrangements within the C-terminal domain of L99A. These rearrangements are identical to those observed along the ground state to excited state transitions characterized by molecular dynamic simulations run on the Anton supercomputer. Analyses of the molecular properties of the residues lining the egress path suggest that protein surface electrostatic potential may play a role in the release mechanism. Simulations of wild-type T4 lysozyme also reveal that benzene-egress-associated dynamics in the L99A mutant are potentially exaggerations of the substrate-processivity-related dynamics of the wild type.

Kinetics of Ligand-Protein Dissociation from All-Atom Simulations: Are We There Yet?

Large parallel gains in the development of both computational resources and sampling methods have now made it possible to simulate dissociation events in ligand-protein complexes with all-atom resolution. Such encouraging progress, together with the inherent spatiotemporal resolution associated with molecular simulations, has left their use for investigating dissociation processes brimming with potential, both in rational drug design, where it can be an invaluable tool for determining the mechanistic driving forces behind dissociation rate constants, and in force-field development, where it can provide a catalog of transient molecular structures with which to refine force fields. Although much progress has been made in making force fields more accurate, reducing their error for transient structures along a transition path could yet prove to be a critical development helping to make kinetic predictions much more accurate. In what follows, we will provide a state-of-the-art compilation of the enhanced sampling methods based on molecular dynamics (MD) simulations used to investigate the kinetics and mechanisms of ligand-protein dissociation processes. Due to the time scales of such processes being slower than what is accessible using straightforward MD simulations, several ingenious schemes are being devised at a rapid rate to overcome this obstacle. Here we provide an up-to-date compendium of such methods and their achievements and shortcomings in extracting mechanistic insight into ligand-protein dissociation. We conclude with a critical and provocative appraisal attempting to answer the title of this Perspective.

Toward Achieving Efficient and Accurate Ligand-Protein Unbinding with Deep Learning and Molecular Dynamics through RAVE

In this work, we demonstrate how to leverage our recent iterative deep learning-all atom molecular dynamics (MD) technique "Reweighted autoencoded variational Bayes for enhanced sampling (RAVE)" (Ribeiro, Bravo, Wang, Tiwary, J. Chem. Phys. 2018, 149, 072301) for investigating ligand-protein unbinding mechanisms and calculating absolute binding free energies, Δ G_b, when plagued with difficult to sample rare events. In order to do so, we introduce a simple but powerful extension to RAVE that allows learning a reaction coordinate expressed as a piecewise function that is linear over all intervals. Such an approach allows us to retain the physical interpretation of a RAVE-derived reaction coordinate while making the method more applicable to a wider range of complex biophysical problems. As we will demonstrate, using as our test-case the slow dissociation of benzene from the L99A variant of lysozyme, the RAVE extension led to observing an unbinding event in 100% of the independent all-atom MD simulations, all within 3-50 ns for a process that takes on an average close to few hundred milliseconds, which reflects a 7 orders of magnitude acceleration relative to straightforward MD. Furthermore, we will show that without the use of time-dependent biasing, clear back-and-forth movement between metastable intermediates was achieved during the various simulations, demonstrating the caliber of the RAVE-derived piecewise reaction coordinate and bias potential, which together drive efficient and accurate sampling of the ligand-protein dissociation event. Last, we report the results for Δ G_b, which via very short MD simulations, can form a strict lower-bound that is ∼2-3 kcal/mol off from experiments. We believe that RAVE, together with its multidimensional extension that we introduce here, will be a useful tool for simulating the slow unbinding process of practical ligand-protein complexes in an automated manner with minimal use of human intuition.

Dynamic disulfide exchange in a crystallin protein in the human eye lens promotes cataract-associated aggregation

Increased light scattering in the eye lens due to aggregation of the long-lived lens proteins, crystallins, is the cause of cataract disease. Several mutations in the gene encoding human γD-crystallin (HγD) cause misfolding and aggregation. Cataract-associated substitutions at Trp⁴² cause the protein to aggregate in vitro from a partially unfolded intermediate locked by an internal disulfide bridge, and proteomic evidence suggests a similar aggregation precursor is involved in age-onset cataract. Surprisingly, WT HγD can promote aggregation of the W42Q variant while itself remaining soluble. Here, a search for a biochemical mechanism for this interaction has revealed a previously unknown oxidoreductase activity in HγD. Using in vitro oxidation, mutational analysis, cysteine labeling, and MS, we have assigned this activity to a redox-active internal disulfide bond that is dynamically exchanged among HγD molecules. The W42Q variant acts as a disulfide sink, reducing oxidized WT and forming a distinct internal disulfide that kinetically traps the aggregation-prone intermediate. Our findings suggest a redox "hot potato" competition among WT and mutant or modified polypeptides wherein variants with the lowest kinetic stability are trapped in aggregation-prone intermediate states upon accepting disulfides from more stable variants. Such reactions may occur in other long-lived proteins that function in oxidizing environments. In these cases, aggregation may be forestalled by inhibiting disulfide flow toward mutant or damaged polypeptides.

Frequency adaptive metadynamics for the calculation of rare-event kinetics

The ability to predict accurate thermodynamic and kinetic properties in biomolecular systems is of both scientific and practical utility. While both remain very difficult, predictions of kinetics are particularly difficult because rates, in contrast to free energies, depend on the route taken. For this reason, specific enhanced sampling methods are needed to calculate long-time scale kinetics. It has recently been demonstrated that it is possible to recover kinetics through the so-called "infrequent metadynamics" simulations, where the simulations are biased in a way that minimally corrupts the dynamics of moving between metastable states. This method, however, requires the bias to be added slowly, thus hampering applications to processes with only modest separations of time scales. Here we present a frequency-adaptive strategy which bridges normal and infrequent metadynamics. We show that this strategy can improve the precision and accuracy of rate calculations at fixed computational cost and should be able to extend rate calculations for much slower kinetic processes.

Kinetics of Huperzine A Dissociation from Acetylcholinesterase via Multiple Unbinding Pathways

The dissociation of huperzine A (hupA) from Torpedo californica acetylcholinesterase ( TcAChE) was investigated by 4 μs unbiased and biased all-atom molecular dynamics (MD) simulations in explicit solvent. We performed our study using memetic sampling (MS) for the determination of reaction pathways (RPs), metadynamics to calculate free energy, and maximum-likelihood estimation (MLE) to recover kinetic rates from unbiased MD simulations. Our simulations suggest that the dissociation of hupA occurs mainly via two RPs: a front door along the axis of the active-site gorge (pwf) and through a new transient side door (pws), i.e., formed by the Ω-loop (residues 67-94 of TcAChE). An analysis of the inhibitor unbinding along the RPs suggests that pws is opened transiently after hupA and the Ω-loop reach a low free-energy transition state characterized by the orientation of the pyridone group of the inhibitor directed toward the Ω-loop plane. Unlike pws, pwf does not require large structural changes in TcAChE to be accessible. The estimated free energies and rates agree well with available experimental data. The dissociation rates along the unbinding pathways are similar, suggesting that the dissociation of hupA along pws is likely to be relevant. This indicates that perturbations to hupA- TcAChE interactions could potentially induce pathway hopping. In summary, our results characterize the slow-onset inhibition of TcAChE by hupA, which may provide the structural and energetic bases for the rational design of the next-generation slow-onset inhibitors with optimized pharmacokinetic properties for the treatment of Alzheimer's disease.

Atomic resolution mechanism of ligand binding to a solvent inaccessible cavity in T4 lysozyme

Ligand binding sites in proteins are often localized to deeply buried cavities, inaccessible to bulk solvent. Yet, in many cases binding of cognate ligands occurs rapidly. An intriguing system is presented by the L99A cavity mutant of T4 Lysozyme (T4L L99A) that rapidly binds benzene (~106 M-1s-1). Although the protein has long served as a model system for protein thermodynamics and crystal structures of both free and benzene-bound T4L L99A are available, the kinetic pathways by which benzene reaches its solvent-inaccessible binding cavity remain elusive. The current work, using extensive molecular dynamics simulation, achieves this by capturing the complete process of spontaneous recognition of benzene by T4L L99A at atomistic resolution. A series of multi-microsecond unbiased molecular dynamics simulation trajectories unequivocally reveal how benzene, starting in bulk solvent, diffuses to the protein and spontaneously reaches the solvent inaccessible cavity of T4L L99A. The simulated and high-resolution X-ray derived bound structures are in excellent agreement. A robust four-state Markov model, developed using cumulative 60 μs trajectories, identifies and quantifies multiple ligand binding pathways with low activation barriers. Interestingly, none of these identified binding pathways required large conformational changes for ligand access to the buried cavity. Rather, these involve transient but crucial opening of a channel to the cavity via subtle displacements in the positions of key helices (helix4/helix6, helix7/helix9) leading to rapid binding. Free energy simulations further elucidate that these channel-opening events would have been unfavorable in wild type T4L. Taken together and via integrating with results from experiments, these simulations provide unprecedented mechanistic insights into the complete ligand recognition process in a buried cavity. By illustrating the power of subtle helix movements in opening up multiple pathways for ligand access, this work offers an alternate view of ligand recognition in a solvent-inaccessible cavity, contrary to the common perception of a single dominant pathway for ligand binding.

Biomolecular conformational changes and ligand binding: from kinetics to thermodynamics

The behaviour of biomolecular systems is governed by their thermodynamic and kinetic properties. It is thus important to be able to calculate, for example, both the affinity and rate of binding and dissociation of a protein-ligand complex, or the populations and exchange rates between distinct conformational states. Because these are typically rare events, calculating these properties from long molecular dynamics simulations remains extremely difficult. Instead, one often adopts a divide-and-conquer strategy in which equilibrium free-energy differences and the fastest state-to-state transition (e.g. ligand association or minor-to-major state conversion) are combined to estimate the slow rate (e.g. ligand dissociation) using a two-state assumption. Here we instead address these problems by using a previously developed method to calculate both the forward and backward rates directly from simulations. We then estimate the thermodynamics from the rates, and validate these values by independent means. We applied the approach to three systems of increasing complexity, including the association and dissociation of benzene to a fully buried cavity inside the L99A mutant variant of T4 lysozyme. In particular, we were able to determine both millisecond association and dissociation rates, and the affinity, of the protein-ligand system by directly observing dozens of rare events in atomic detail. Our approach both sheds light on the precision of methods for calculating kinetics and further provides a generally useful test for the internal consistency of kinetics and thermodynamics. We also expect our route to be useful for obtaining both the kinetics and thermodynamics at the same time in more challenging cases.

Characterizing Drug-Target Residence Time with Metadynamics: How To Achieve Dissociation Rate Efficiently without Losing Accuracy against Time-Consuming Approaches

Drug-target residence time plays a vital role in drug efficacy. However, there is still no effective strategy to predict drug residence time. Here, we propose to use the optimized (or minimized) structures derived from holo-state proteins to calculate drug residence time, which could give a comparable or even better prediction accuracy compared with those calculated utilizing a large number of molecular dynamics (MD) structures based on the Poisson process. Besides, in addition to the Poisson process, one may use fewer samples for predicting residence time due to the reason that, in a large extent, the calculated drug residence time is stable and independent of the number of samples used for the prediction. With remarkably reduced computational load, the proposed strategy may be promising for large-scale drug residence time prediction, such as post-processing in virtual screening (VS) and lead compound optimization.

An optimal distance cutoff for contact-based Protein Structure Networks using side-chain centers of mass

Proteins are highly dynamic entities attaining a myriad of different conformations. Protein side chains change their states during dynamics, causing clashes that are propagated at distal sites. A convenient formalism to analyze protein dynamics is based on network theory using Protein Structure Networks (PSNs). Despite their broad applicability, few efforts have been devoted to benchmarking PSN methods and to provide the community with best practices. In many applications, it is convenient to use the centers of mass of the side chains as nodes. It becomes thus critical to evaluate the minimal distance cutoff between the centers of mass which will provide stable network properties. Moreover, when the PSN is derived from a structural ensemble collected by molecular dynamics (MD), the impact of the MD force field has to be evaluated. We selected a dataset of proteins with different fold and size and assessed the two fundamental properties of the PSN, i.e. hubs and connected components. We identified an optimal cutoff of 5 Å that is robust to changes in the force field and the proteins. Our study builds solid foundations for the harmonization and standardization of the PSN approach.

How and when does an anticancer drug leave its binding site?

Obtaining atomistic resolution of drug unbinding from a protein is a much sought-after experimental and computational challenge. We report the unbinding dynamics of the anticancer drug dasatinib from c-Src kinase in full atomistic resolution using enhanced sampling molecular dynamics simulations. We obtain multiple unbinding trajectories and determine a residence time in agreement with experiments. We observe coupled protein-water movement through multiple metastable intermediates. The water molecules form a hydrogen bond bridge, elongating a specific, evolutionarily preserved salt bridge and enabling conformation changes essential to ligand unbinding. This water insertion in the salt bridge acts as a molecular switch that controls unbinding. Our findings provide a mechanistic rationale for why it might be difficult to engineer drugs targeting certain specific c-Src kinase conformations to have longer residence times.

Ligand diffusion in proteins via enhanced sampling in molecular dynamics

Computational simulations in biophysics describe the dynamics and functions of biological macromolecules at the atomic level. Among motions particularly important for life are the transport processes in heterogeneous media. The process of ligand diffusion inside proteins is an example of a complex rare event that can be modeled using molecular dynamics simulations. The study of physical interactions between a ligand and its biological target is of paramount importance for the design of novel drugs and enzymes. Unfortunately, the process of ligand diffusion is difficult to study experimentally. The need for identifying the ligand egress pathways and understanding how ligands migrate through protein tunnels has spurred the development of several methodological approaches to this problem. The complex topology of protein channels and the transient nature of the ligand passage pose difficulties in the modeling of the ligand entry/escape pathways by canonical molecular dynamics simulations. In this review, we report a methodology involving a reconstruction of the ligand diffusion reaction coordinates and the free-energy profiles along these reaction coordinates using enhanced sampling of conformational space. We illustrate the above methods on several ligand-protein systems, including cytochromes and G-protein-coupled receptors. The methods are general and may be adopted to other transport processes in living matter.