Dickson Lab @ MSU

Targeted protein degradation is an emerging approach that utilizes cellular degradation pathways to inhibit a target protein. Small molecules such as molecular glues or PROTACs can be used to mediate the formation of a ternary complex with an E3 ligase and the target protein, which can dramatically enhance the degradation process. This approach is promising for cancer therapy, where degradation of oncogenic proteins can lead to cancer cell toxicity. To design new molecular glues, it is important to develop methods that predict how well a given molecule stabilizes a protein-protein interaction. However, conventional molecular dynamics simulations...

Neurolysin (Nln) is an M3 metallopeptidase that regulates neuropeptide concentration in the central nervous system. It has emerged as a therapeutic target for mitigating post-ischemic injury by hydrolyzing and inactivating several neuropeptides. It has been recently shown that small molecule activators, such as pyridine piperazine (Py-Pip) derivatives, can enhance Nln catalytic activity, facilitating hydrolysis of Nln substrate peptides. However, binding sites of these molecules and the mechanism of action remain unclear due to the dynamic nature of Nln. Here, we use molecular dynamics (MD) simulation of the apo and activator-bound Nln systems along with Markov state...

Peptides are desirable therapeutics due to their inherent potency, safety, cost-effectiveness and ability to engage large or more complex protein surfaces. Slower kinetics of protein-peptide (un)binding can directly influence their drug efficacy and duration of action, in part by improving plasma stability of the peptide. The CLR:RAMP1 complex and its endogenous agonist peptide CGRP are of particularly high interest because of their central role in migraine pathophysiology. A better understanding of peptide binding mechanisms is needed for the development of next-generation peptide-based drugs with optimized kinetic properties. In this study, we comparatively analyze constructs of native...

A common first step in drug design is virtual high throughput screening (VHTS), where a large number of potential drug molecules are computationally modeled in a protein binding pocket and filtered down to a smaller set of hits that can be further tested computationally or experimentally. Traditional strategies for VHTS do not account for ligand-induced conformational changes in proteins, as they typically rely on a single static structure to represent the protein. This neglects the role of binding entropy and the fact that different ligand molecules can induce slightly different conformations in the protein binding site...

Calcitonin gene-related peptide (CGRP) is a 37-amino acid neuropeptide that functions in pain signaling and neuroimmune communication. The CGRP receptor, CGRPR, is a class B GPCR that is a drug target for migraine headache and other disorders. Here, we used nanoBRET receptor binding and cAMP biosensor signaling assays and theoretical modeling to characterize the CGRPR “two-domain” peptide binding mechanism. Single-site extracellular domain (ECD)-binding and two-site ECD/transmembrane domain (TMD)-binding peptides were examined for CGRP and a high-affinity variant “ssCGRP” with modifications in the C-terminal region. Wildtype and ssCGRP(27-37) bound the ECD with affinities of 1 μM and...

Accurate prediction of molecular geometries is crucial for drug discovery and materials science. Existing fast conformer prediction algorithms often rely on approximate empirical energy functions, resulting in low accuracy. More accurate methods like ab initio molecular dynamics and Markov chain Monte Carlo can be computationally expensive due to the need for evaluating quantum mechanical energy functions. To address this, we introduce AGDIFF, a novel machine learning framework that utilizes diffusion models for efficient and accurate molecular structure prediction. AGDIFF extends previous models (such as GeoDiff) by enhancing the global, local, and edge encoders with attention mechanisms,...

The weighted ensemble (WE) algorithm is gaining popularity as a rare event method for studying long timescale processes with molecular dynamics. WE is particularly useful for determining kinetic properties, like rates of protein (un)folding and ligand (un)binding, where transition rates can be calculated from the flux of trajectories into a target basin of interest. However, this flux depends exponentially on the number of splitting events that a given trajectory experiences before reaching the target state, and can vary by orders of magnitude between WE replicates. Markov state models (MSM) are helpful tools to aggregate information across...

For many drug targets, it has been shown that the kinetics of drug binding (e.g., on rate and off rate) is more predictive of drug efficacy than thermodynamic quantities alone. This motivates the development of predictive computational models that can be used to optimize compounds on the basis of their kinetics. The structural details underpinning these computational models are found not only in the bound state but also in the short-lived ligand binding transition states. Although transition states cannot be directly observed experimentally due to their extremely short lifetimes, recent successes have demonstrated that modeling the...

Ligand design problems involve searching chemical space for a molecule with a set of desired properties. As chemical space is discrete, this search must be conducted in a pointwise manner, separately investigating one molecule at a time, which can be inefficient. We propose a method called “Flexible Topology”, where a ligand is composed of a set of shapeshifting “ghost” atoms, whose atomic identities and connectivity can dynamically change over the course of a simulation. Ghost atoms are guided toward their target positions using a translation-, rotation-, and index-invariant restraint potential. This is the first step toward...

The signaling peptides adrenomedullin 2/intermedin (AM2/IMD), adrenomedullin (AM), and CGRP have overlapping and distinct functions in the cardiovascular, lymphatic, and nervous systems by activating three shared receptors comprised of the class B GPCR CLR in complex with a RAMP1, −2, or −3 modulatory subunit. Here, we report that AM2/IMD, which is thought to be a non-selective agonist, is kinetically selective for CLR-RAMP3, known as the AM2R. AM2/IMD-AM2R elicited substantially longer duration cAMP signaling than the eight other peptide-receptor combinations due to AM2/IMD slow off-rate binding kinetics. The regions responsible for the slow off-rate were mapped to...

The prediction of (un)binding rates and free energies is of great significance to the drug design process. Although many enhanced sampling algorithms and approaches have been developed, there is not yet a reliable workflow to predict these quantities. Previously we have shown that free energies and transition rates can be calculated by directly simulating the binding and unbinding processes with our variant of the WE algorithm “Resampling of Ensembles by Variation Optimization”, or “REVO”. Here, we calculate binding free energies retrospectively for three SAMPL6 host-guest systems and prospectively for a SAMPL9 system to test a modification...

Targeted protein degradation (TPD) is a promising approach in drug discovery for degrading proteins implicated in diseases. A key step in this process is the formation of a ternary complex where a heterobifunctional molecule induces proximity of an E3 ligase to a protein of interest (POI), thus facilitating ubiquitin transfer to the POI. In this work, we characterize 3 steps in the TPD process. (1) We simulate the ternary complex formation of SMARCA2 bromodomain and VHL E3 ligase by combining hydrogen-deuterium exchange mass spectrometry with weighted ensemble molecular dynamics (MD). (2) We characterize the conformational heterogeneity...

Chaperones and other quality control machinery guard proteins from inappropriate aggregation, which is a hallmark of neurodegenerative diseases. However, how the systems that regulate the foldedness of the proteome remain buffered under stress conditions and in different cellular compartments remains incompletely understood. In this study, we applied a FRET-based strategy to explore how well quality control machinery protects against the misfolding and aggregation of bait biosensor proteins, made from the prokaryotic ribonuclease barnase, in the nucleus and cytosol of human embryonic kidney 293T cells. We found that those barnase biosensors were prone to misfolding, were less...

Improper reaction coordinates can pose significant problems for path-based binding free energy calculations. Particularly, omission of long timescale motions can lead to over-estimation of the energetic barriers between the bound and unbound states. Many methods exist to construct the optimal reaction coordinate using a pre-defined basis set of features. Although simulations are typically conducted in explicit solvent, the solvent atoms are often excluded by these feature sets—resulting in little being known about their role in reaction coordinates, and ultimately, their role in determining (un)binding rates and free energies. In this work, analysis is done on an...

The human ACE2 enzyme serves as a critical first recognition point of coronaviruses, including SARS-CoV-2. In particular, the extracellular domain of ACE2 interacts directly with the S1 tailspike protein of the SARS-CoV-2 virion through a broad protein–protein interface. Although this interaction has been characterized by X-ray crystallography, these structures do not reveal significant differences in the ACE2 structure upon S1 protein binding. In this work, using several all-atom molecular dynamics simulations, we show persistent differences in the ACE2 structure upon binding. These differences are determined with the linear discriminant analysis (LDA) machine learning method and validated...

The prediction of log P values is one part of the statistical assessment of the modeling of proteins and ligands (SAMPL) blind challenges. Here, we use a molecular graph representation method called Geometric Scattering for Graphs (GSG) to transform atomic attributes to molecular features. The atomic attributes used here are parameters from classical molecular force fields including partial charges and Lennard–Jones interaction parameters. The molecular features from GSG are used as inputs to neural networks that are trained using a “master” dataset comprised of over 41,000 unique log P values. The specific molecular targets in the...

This work examines methods for predicting the partition coefficient (logP) for a dataset of small molecules. Here, we use atomic attributes such as radius and partial charge, which are typically used as force field parameters in classical molecular dynamics simulations. These atomic features are transformed into index-invariant molecular features using a recently developed method called Geometric Scattering for Graphs (GSG). We call this approach “ClassicalGSG” and examine its performance under a broad range of conditions and hyperparameters. We train a ClassicalGSG logP predictor with neural networks using 10,722 molecules from the ChEMBL21 dataset and apply it...

The translocator protein (TSPO), previously known as the peripheral benzodiazepine receptor, is of longstanding medical interest as both a biomarker for neuroinjury and a potential drug target for neuroinflammation and other disorders. Recently it was shown that ligand residence time is a key factor determining steroidogenic efficacy of TSPO-binding compounds. This spurs interest in simulations of (un)binding pathways of TSPO ligands, which could reveal the molecular interactions governing ligand residence time. In this study, we use a weighted ensemble algorithm to determine the unbinding pathway for different poses of PK-11195, a TSPO ligand used in neuroimaging....

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...

The free energy of a process is the fundamental quantity that determines its spontaneity or propensity at a given temperature. In particular, the binding free energy of a drug candidate to its biomolecular target is used as an objective quantity in drug design. Recently, binding kinetics—rates of association (kon) and dissociation (koff)—have also demonstrated utility for their ability to predict efficacy and in some cases have been shown to be more predictive than the binding free energy alone. Some methods exist to calculate binding kinetics from molecular simulations, although these are typically more difficult to calculate...

Approaches for computing small molecule binding free energies based on molecular simulations are now regularly being employed by academic and industry practitioners to study receptor-ligand systems and prioritize the synthesis of small molecules for ligand design. Given the variety of methods and implementations available, it is natural to ask how the convergence rates and final predictions of these methods compare. In this study, we describe the concept and results for the SAMPL6 SAMPLing challenge, the first challenge from the SAMPL series focusing on the assessment of convergence properties and reproducibility of binding free energy methodologies. We...

The free energy of transitions between stable states is the key thermodynamic quantity that governs the relative probabilities of the forward and reverse reactions and the ratio of state probabilities at equilibrium. The binding free energy of a drug and its receptor is of particular interest, as it serves as an optimization function for drug design. Over the years, many computational methods have been developed to calculate binding free energies, and while many of these methods have a long history, issues such as convergence of free energy estimates and the projection of a binding process onto...

The weighted ensemble (WE) strategy has been demonstrated to be highly efficient in generating pathways and rate constants for rare events such as protein folding and protein binding using atomistic molecular dynamics simulations. Here we present five tutorials instructing users in the best practices for preparing, carrying out, and analyzing WE simulations for various applications using the WESTPA software. Users are expected to already have significant experience with running standard molecular dynamics simulations using the underlying dynamics engine of interest (e.g. Amber, Gromacs, OpenMM). The tutorials range from a molecular association process in explicit solvent to...

Conventional molecular dynamics simulations are incapable of sampling many important interactions in biomolecular systems due to their high dimensionality and rough energy landscapes. To observe rare events and calculate transition rates in these systems, enhanced sampling is a necessity. In particular, the study of ligand-protein interactions necessitates a diverse ensemble of protein conformations and transition states, and for many systems, this occurs on prohibitively long time scales. Previous strategies such as WExplore that can be used to determine these types of ensembles are hindered by problems related to the regioning of conformational space. Here, we propose...

Electrostatic channeling is a naturally-occurring approach to control the flux of charged intermediates in catalytic cascades. Computational techniques have enabled quantitative understanding of such mechanisms, augmenting experimental approaches by modeling molecular interactions in atomic detail. In this work, we report the first utilization of a Markov State Model (MSM) to describe the surface diffusion of a reaction intermediate, glucose 6-phosphate, on an artificially modified cascade where hexokinase and glucose-6-phosphate dehydrogenase are covalently conjugated by a cationic oligopeptide bridge. Conformation space networks are used to represent intermediate transport on enzyme surfaces, along with committor probabilities that assess...

Bromodomain and PHD finger containing protein transcription factor (BPTF) is an epigenetic protein involved in chromatin remodelling and is a potential anticancer target. The BPTF bromodomain has one reported small molecule inhibitor (AU1, rac-1). Here, advances made on the study of rac-1 are reported, specifically, identification of the active enantiomer for binding to the bromodomain, and reduction in off-target binding to known kinase targets. Additionally, a ligand deconstruction analysis was conducted to characterize important pharmacophores for engaging the BPTF bromodomain. These studies have been enabled by a protein-based fluorine NMR approach, highlighting the versatility of the...

The interaction between a ligand and a protein involves a multitude of conformational states. To achieve a particular deeply-bound pose the ligand must search across a rough free energy landscape, with many metastable minima. Creating maps of the ligand binding landscape is a great challenge, as binding and release events typically occur on timescales that are beyond the reach of molecular simulation. The WExplore enhanced sampling method is well-suited to build these maps, as it is designed to broadly explore free-energy landscapes, and is capable of simulating ligand release pathways that occur on timescales as long...

Interest in ligand binding kinetics has been growing rapidly, as it is being discovered in more and more systems that ligand residence time is the crucial factor governing drug efficacy. Many enhanced sampling methods have been developed with the goal of predicting ligand binding rates (kon) and/or ligand unbinding rates (koff) through explicit simulation of ligand binding pathways, and these methods work by very different mechanisms. Although there is not yet a blind challenge for ligand binding kinetics, here we take advantage of experimental measurements and rigorously computed benchmarks to compare estimates of KD calculated as...

Signaling molecule phosphatidylinositol 4,5-bisphosphate (PIP2) is produced primarily by phosphatidylinositol 4-phosphate 5-kinase (PIP5K). PIP5K is essential for the development of the human neuronal system, which has been exemplified by a recessive genetic disorder, lethal congenital contractural syndrome type 3 (LCCS3), caused by a single aspartate-to-asparagine mutation in the kinase domain of PIP5Ky. So far, the exact role of this aspartate residue has yet to be elucidated. In this work, we conducted structural, functional and computational studies on a zebrafish PIP5Ka variant with a mutation at the same site. Compared with the structure of the wild type...

The pool of quality control proteins (QC) that maintains protein-folding homeostasis (proteostasis) is dynamic but can become depleted in human disease. A challenge has been in quantitatively defining the depth of the QC pool. With a new biosensor, flow cytometry-based methods and mathematical modeling we measure the QC capacity to act as holdases and suppress biosensor aggregation. The biosensor system comprises a series of barnase kernels with differing folding stability that engage primarily with HSP70 and HSP90 family proteins. Conditions of proteostasis stimulation and stress alter QC holdase activity and aggregation rates. The method reveals the...

Ligand (un)binding kinetics is being recognized as a determinant of drug specificity and efficacy in an increasing number of systems. However, the calculation of kinetics and the simulation of drug unbinding is more difficult than computing thermodynamic quantities, such as binding free energies. Here we present the first full simulations of an unbinding process at pharmacologically relevant timescales (11 min), without the use of biasing forces, detailed prior knowledge, or specialized processors using the weighted ensemble based algorithm, WExplore. These simulations show the inhibitor TPPU unbinding from its enzyme target soluble epoxide hydrolase (sEH), which is...

Crystal structures of neurolysin, which is a zinc metallopeptidase (neuropeptidase), do not show significant conformational changes upon the binding of an allosteric inhibitor. Neurolysin has a prolate ellipsoid shape with a deep channel that runs almost the entire length of the molecule. In this deep channel, neurolysin hydrolyzes the short neuropeptide neurotensin to create inactive shorter fragments and thus controls the neurotensin level in the tissue. The protein is of interest as a therapeutic target since changes in neurotensin level have been implicated in cardiovascular and neurological disorders and cancer, and inhibitors of neurolysin activity have...

Ligand residence times and binding rates have been found to be useful quantities to consider during drug design. The underlying structural and dynamic determinants of these kinetic quantities are difficult to discern. Driven by developments in computational hardware and simulation methodologies, molecular dynamics (MD) studies of full binding and unbinding pathways have emerged recently, showing these structural and dynamic determinants in atomic detail. However, the long timescales related to drug binding and release are still prohibitive to conventional MD simulation. Here we discuss a suite of enhanced sampling methods that have been applied to the study...

We report simulations of full ligand exit pathways for the trypsin-benzamidine system, generated using the sampling technique WExplore. WExplore is able to observe millisecond-scale unbinding events using many nanosecond-scale trajectories that are run without introducing biasing forces. The algorithm generates rare events by dividing the coordinate space into regions, on-the-fly, and balancing computational effort between regions through cloning and merging steps, as in the weighted ensemble method. The averaged exit flux yields a ligand exit rate of 180 microseconds, which is within an order of magnitude of the experimental value. We obtain broad sampling of ligand...

Proteins can be destabilized by a number of environmental factors such as temperature, pH and mutation. The ability to restore function by small molecule stabilizers, or the introduction of disulde bonds, would be a very powerful tool, but the physical principles that drive this stabilization are not well understood. The first problem lies is in choosing an appropriate binding site or disulfide bond location that will best confer stability to the active site and restore function. Here we present a general framework for predicting which allosteric binding sites correlate with stability in the active site. Using...

Chaperones maintain a healthy proteome by preventing aggregation and by aiding in protein folding. Precisely how chaperones influence the conformational properties of their substrates, however, remains unclear. To achieve a detailed description of dynamic chaperone-substrate interactions, we fused site-specific NMR information with coarse-grained simulations. Our model system is the binding and folding of a chaperone substrate, immunity protein 7 (Im7), with the chaperone Spy. We first used an automated procedure in which NMR chemical shifts inform the construction of system-specific force fields that describe each partner individually. The models of the two binding partners are then...

The binding of ligands with their molecular receptors is of tremendous importance in biology. Although much emphasis has been placed on characterizing binding sites and bound poses that determine the binding thermodynamics, the pathway by which a ligand binds importantly determines the binding kinetics. The computational study of entire unbiased ligand binding and release pathways is still an emerging field, made possible only recently by advances in computational hardware and sampling methodologies. We have developed one such method (WExplore) that is based on a weighted ensemble of trajectories, which we apply to ligand release for the...

Intrinsically disordered regions of proteins can gain structure by binding to a partner. This process, of coupled folding and binding (CFaB), is a fundamental part of many important biological processes. Structure-based models have proven themselves capable of revealing fundamental aspects of how CFaB occurs, however, typical methods to enhance the sampling of these transitions, such as replica exchange, do not adequately sample the transition state region of this extremely rare process. Here, we use a variant of Umbrella Sampling to enforce sampling of the transition states of CFaB of HdeA monomers at neutral pH, an extremely...

The pH-sensing chaperone HdeA promotes the survival of enteropathogenic bacteria during transit through the harshly acidic environment of the mammalian stomach. At low pH, HdeA transitions from an inactive, folded, dimer to chaperone-active, disordered, monomers to protect against the acidiinduced aggregation of periplasmic proteins. Toward achieving a detailed mechanistic understanding of the pH response of HdeA, we develop a multiscale modeling approach to capture its pH-dependent thermodynamics. Our approach combines pKa (logarithmic acid dissociation constant) calculations from all-atom constant pH molecular dynamics simulations with coarse-grained modeling and yields new, atomic-level, insights into HdeA chaperone function that...

A recently engineered mutant of cyan fluorescent protein (WasCFP) that exhibits pH-dependent absorption suggests that its tryptophan-based chromophore switches between neutral (protonated) and charged (deprotonated) states depending on external pH. At pH 8.1, the latter gives rise to green fluorescence as opposed to the cyan color of emission that is characteristic for the neutral form at low pH. Given the high energy cost of deprotonating the tryptophan at the indole nitrogen, this behavior is puzzling, even if the stabilizing effect of the V61K mutation in proximity to the protonation/deprotonation site is considered. Because of its potential...

HIV-1 TAR RNA is a two-helix bulge motif that plays a critical role in HIV viral replication and is an important drug target. However, efforts at designing TAR inhibitors have been challenged by its high degree of structural flexibility, which includes slow large-amplitude reorientations of its helices with respect to one another. Here, we use the recently introduced algorithm WExplore in combination with Euler angles to achieve unprecedented sampling of the TAR conformational ensemble. Our ensemble achieves similar agreement with experimental NMR data when compared with previous TAR computational studies, and is generated at a fraction...

As most relevant motions in biomolecular systems are inaccessible to conventional molecular dynamics simulations, algorithms that enhance sampling of rare events are indispensable. Increasing interest in intrinsically disordered systems and the desire to target ensembles of protein conformations (rather than single structures) in drug development motivate the need for enhanced sampling algorithms that are not limited to two-basin problems, and can efficiently determine structural ensembles. For systems that are not well-studied, this must often be done with little or no information about the dynamics of interest. Here we present a novel strategy to determine structural ensembles...

For cells to function, the concentrations of all proteins in the cell must be maintained at the proper levels (proteostasis). This task - complicated by cellular stresses, protein misfolding, aggregation, and degradation - is performed by a collection of chaperones that alter the configurational landscape of a given client protein through the formation of protein-chaperone complexes. The set of all such complexes and the transitions between them form the proteostasis network. Recently, a computational model was introduced (FoldEco) that synthesizes experimental data into a system-wide description of the proteostasis network of E. coli. This model describes...

The small pH stress-sensing chaperone HdeA helps pathogenic enteric E. coli survive passage through the severely acidic environment of the mammalian stomach. Under stress conditions, HdeA transitions from an inactive folded dimer to a chaperone-active unfolded monomer to prevent the acid-induced aggregation of periplasmic proteins. Here we use a topology-based Go-like model to delineate the relationship between dimer interface formation and monomer folding and to better understand the structural details of the chaperone activation mechanism. Free energy surfaces show that dimer interface formation and monomer folding proceed concurrently through an on-pathway dimeric intermediate in which one...

It has been suggested that the native state of a protein acts as a kinetic hub that can facilitate transitions between nonnative states. Using recently developed tools to quantify mediation probabilities (hub scores), we quantify hub-like behavior in atomic resolution trajectories for the first time. We use a data set of trajectory ensembles for 12 fast-folding proteins previously published by D. E. Shaw Research (Lindorff-Larsen, K.; et al. How Fast-Folding Proteins Fold. Science 2011, 334, 517) with an aggregate simulation time of over 8.2 ms. We visualize the free-energy landscape of each molecule using configuration space...

The free energy landscape of a protein is a function of many interdependent degrees of freedom. For this reason, conceptual constructs (e.g., funnels) have been useful to visualize these landscapes. One relatively new construct is the idea of a hub-like native state that is the final destination of many noninterconverting folding pathways. This is in contrast to the idea of a single predominant folding pathway connecting the native state to a rapidly interconverting ensemble of unfolded states. The key quantity to distinguish between these two ideas is the connectivity of the unfolded ensemble. We present a...

Large deviation theory has emerged as the natural language with which to study transitions between dynamical phases, such as the glass transition. Here we use this approach to show that the entrainment of an oscillator to an external driving force can be described as a phase transition between two regions in a joint space-time representation. Specifically, we numerically obtain exact solutions of the large deviation function for a discrete, finite model of an oscillator under a periodically varying external force. We find that the first derivative of the expectation value of the current diverges in the...