Adaptive BP-Dock: An Induced Fit Docking Approach for Full Receptor Flexibility

DSDPFlex: Flexible-Receptor Docking with GPU Acceleration

Molecular docking is an essential tool in structure-based drug discovery, widely utilized to model ligand-protein interactions and enrich potential hits. Among the different docking strategies, semiflexible docking (rigid-receptor and flexible-ligand model) is the most popular, benefiting from its balance of docking accuracy and speed. However, this approach ignores the conformational changes of proteins and hence demands suitable protein conformations as input. When the binding interaction adheres to an induced-fit model, flexible methods such as molecular dynamics simulation can be utilized, but they are computationally demanding. To balance between speed and accuracy, the flexible docking approach is an effective choice, as exemplified by AutoDock Vina and AutoDockFR, which treat selected protein side chains as flexible parts. However, the efficiency of flexible docking methods is yet to be improved for virtual screening usage. In this article, we introduce DSDPFlex, an improved flexible-receptor docking method accelerated by GPU parallelization. Beyond acceleration, optimizations with respect to sampling, scoring, and search space are implemented in DSDPFlex to further improve its capability in flexible tasks. In cross-docking evaluation, DSDPFlex demonstrates superior accuracy compared to AutoDock Vina and is 100 times faster than Vina in flexible-receptor tasks. We also show the advantage of flexible-receptor methods on suboptimal pockets and validate the advantage of DSDPFlex in screening on apo and AlphaFold2-predicted structures. With improvements in both efficiency and accuracy, DSDPFlex is expected to hold potential in future docking-based studies.

Engineering gain-of-function mutants of a WW domain by dynamics and structural analysis

Proteins gain optimal fitness such as foldability and function through evolutionary selection. However, classical studies have found that evolutionarily designed protein sequences alone cannot guarantee foldability, or at least not without considering local contacts associated with the initial folding steps. We previously showed that foldability and function can be restored by removing frustration in the folding energy landscape of a model WW domain protein, CC16, which was designed based on Statistical Coupling Analysis (SCA). Substitutions ensuring the formation of five local contacts identified as "on-path" were selected using the closest homolog native folded sequence, N21. Surprisingly, the resulting sequence, CC16-N21, bound to Group I peptides, while N21 did not. Here, we identified single-point mutations that enable N21 to bind a Group I peptide ligand through structure and dynamic-based computational design. Comparison of the docked position of the CC16-N21/ligand complex with the N21 structure showed that residues at positions 9 and 19 are important for peptide binding, whereas the dynamic profiles identified position 10 as allosterically coupled to the binding site and exhibiting different dynamics between N21 and CC16-N21. We found that swapping these positions in N21 with matched residues from CC16-N21 recovers nature-like binding affinity to N21. This study validates the use of dynamic profiles as guiding principles for affecting the binding affinity of small proteins.

Computation of the Protein Conformational Transition Pathway on Ligand Binding by Linear Response-Driven Molecular Dynamics

While extremely important for relating the protein structure to its biological function, determination of the protein conformational transition pathway upon ligand binding is made difficult due to the transient nature of intermediates, a large and rugged conformational space, and coupling between protein dynamics and ligand-protein interactions. Existing methods that rely on prior knowledge of the bound (holo) state structure are restrictive. A second concern relates to the correspondence of intermediates obtained to the metastable states on the apo → holo transition pathway. Here, we have taken the protein apo structure and ligand-binding site as only inputs and combined an elastic network model (ENM) representation of the protein Hamiltonian with linear response theory (LRT) for protein-ligand interactions to identify the set of slow normal modes of protein vibrations that have a high overlap with the direction of the protein conformational change. The structural displacement along the chosen direction was performed using excited normal modes molecular dynamics (MDeNM) simulations rather than by the direct use of LRT. Herein, the MDeNM excitation velocity was optimized on-the-fly on the basis of its coupling to protein dynamics and ligand-protein interactions. Thus, a determined set of structures was validated against crystallographic and simulation data on four protein-ligand systems, namely, adenylate kinase-di(adenosine-5')pentaphosphate, ribose binding protein-β-d-ribopyranose, DNA β-glucosyltransferase-uridine-5'-diphosphate, and G-protein α subunit-guanosine-5'-triphosphate, which present important differences in protein conformational heterogeneity, ligand binding mechanism, viz. induced-fit or conformational selection, extent, and nonlinearity in protein conformational changes upon ligand binding, and presence of allosteric effects. The obtained set of intermediates was used as an input to path metadynamics simulations to obtain the free energy profile for the apo → holo transition.

CIFDock: A novel CHARMM-based flexible receptor-flexible ligand docking protocol

Docking studies play a critical role in the current workflow of drug discovery. However, limitations may often arise through factors including inadequate ligand sampling, a lack of protein flexibility, scoring function inadequacies (e.g., due to metals, co-factors, etc.), and difficulty in retaining explicit water molecules. Herein, we present a novel CHARMM-based induced fit docking (CIFDock) workflow that can circumvent these limitations by employing all-atom force fields coupled to enhanced sampling molecular dynamics procedures. Self-guided Langevin dynamics simulations are used to effectively sample relevant ligand conformations, side chain orientations, crystal water positions, and active site residue motion. Protein flexibility is further enhanced by dynamic sampling of side chain orientations using an expandable rotamer library. Steps in the procedure consisting of fixing individual components (e.g., the ligand) while sampling the other components (e.g., the residues in the active site of the protein) allow for the complex to adapt to conformational changes. Ultimately, all components of the complex-the protein, ligand, and waters-are sampled simultaneously and unrestrained with SGLD to capture any induced fit effects. This modular flexible docking procedure is automated using CHARMM scripting, interfaced with SLURM array processing, and parallelized to use the desired number of processors. We validated the CIFDock procedure by performing cross-docking studies using a data set comprised of 21 pharmaceutically relevant proteins. Five variants of the CHARMM-based SWISSDOCK scoring functions were created to quantify the results of the final generated poses. Results obtained were comparable to, or in some cases improved upon, commercial docking program data.

Design of novel cyanovirin-N variants by modulation of binding dynamics through distal mutations

We develop integrated co-evolution and dynamic coupling (ICDC) approach to identify, mutate, and assess distal sites to modulate function. We validate the approach first by analyzing the existing mutational fitness data of TEM-1 β-lactamase and show that allosteric positions co-evolved and dynamically coupled with the active site significantly modulate function. We further apply ICDC approach to identify positions and their mutations that can modulate binding affinity in a lectin, cyanovirin-N (CV-N), that selectively binds to dimannose, and predict binding energies of its variants through Adaptive BP-Dock. Computational and experimental analyses reveal that binding enhancing mutants identified by ICDC impact the dynamics of the binding pocket, and show that rigidification of the binding residues compensates for the entropic cost of binding. This work suggests a mechanism by which distal mutations modulate function through dynamic allostery and provides a blueprint to identify candidates for mutagenesis in order to optimize protein function.

Preliminary modelling as the first stage of targeted organic synthesis

The review aims to present a classification and applicability analysis of methods for preliminary molecular modelling for targeted organic, catalytic and biocatalytic synthesis. The following three main approaches are considered as a primary classification of the methods: modelling of the target – ligand coordination without structural information on both the target and the resulting complex; calculations based on experimentally obtained structural information about the target; and dynamic simulation of the target – ligand complex and the reaction mechanism with calculation of the free energy of the reaction. The review is meant for synthetic chemists to be used as a guide for building an algorithm for preliminary modelling and synthesis of structures with specified properties.

The bibliography includes 353 references.

Substitutions at Nonconserved Rheostat Positions Modulate Function by Rewiring Long-Range, Dynamic Interactions

Amino acid substitutions at nonconserved protein positions can have noncanonical and "long-distance" outcomes on protein function. Such outcomes might arise from changes in the internal protein communication network, which is often accompanied by changes in structural flexibility. To test this, we calculated flexibilities and dynamic coupling for positions in the linker region of the lactose repressor protein. This region contains nonconserved positions for which substitutions alter DNA-binding affinity. We first chose to study 11 substitutions at position 52. In computations, substitutions showed long-range effects on flexibilities of DNA-binding positions, and the degree of flexibility change correlated with experimentally measured changes in DNA binding. Substitutions also altered dynamic coupling to DNA-binding positions in a manner that captured other experimentally determined functional changes. Next, we broadened calculations to consider the dynamic coupling between 17 linker positions and the DNA-binding domain. Experimentally, these linker positions exhibited a wide range of substitution outcomes: Four conserved positions tolerated hardly any substitutions ("toggle"), ten nonconserved positions showed progressive changes from a range of substitutions ("rheostat"), and three nonconserved positions tolerated almost all substitutions ("neutral"). In computations with wild-type lactose repressor protein, the dynamic couplings between the DNA-binding domain and these linker positions showed varied degrees of asymmetry that correlated with the observed toggle/rheostat/neutral substitution outcomes. Thus, we propose that long-range and noncanonical substitutions outcomes at nonconserved positions arise from rewiring long-range communication among functionally important positions. Such calculations might enable predictions for substitution outcomes at a range of nonconserved positions.

Atlasing of Assembly Landscapes using Distance Geometry and Graph Rigidity

This Article describes a novel geometric methodology for analyzing free energy and kinetics of assembly driven by short-range pair-potentials in an implicit solvent and provides a proof-of-concept illustration of its unique capabilities. An atlas is a labeled partition of the assembly landscape into a roadmap of maximal, contiguous, nearly-equipotential-energy conformational regions or macrostates, together with their neighborhood relationships. The new methodology decouples the roadmap generation from sampling and produces: (1) a queryable atlas of local potential energy minima, their basin structure, energy barriers, and neighboring basins; (2) paths between a specified pair of basins, each path being a sequence of conformational regions or macrostates below a desired energy threshold; and (3) approximations of relative path lengths, basin volumes (configurational entropy), and path probabilities. Results demonstrating the core algorithm's capabilities and high computational efficiency have been generated by a resource-light, curated open source software implementation EASAL (Efficient Atlasing and Search of Assembly Landscapes, ACM Trans. Math. Softw. 2018 44, 1-48. 10.1145/3204472; see software, Efficient Atlasing and Search of Assembly Landscapes, 2016. https://bitbucket.org/geoplexity/easal; video, Video Illustrating the opensource software EASAL, 2016. https://cise.ufl.edu/~sitharam/EASALvideo.mpeg; and user guide, EASAL software user guide, 2016. https://bitbucket.org/geoplexity/easal/src/master/CompleteUserGuide.pdf). Running on a laptop with Intel(R) Core(TM) i7-7700@3.60 GHz CPU with 16GB of RAM, EASAL atlases several hundred thousand conformational regions or macrostates in minutes using a single compute core. Subsequent path and basin computations each take seconds. A parallelized EASAL version running on the same laptop with 4 cores gives a 3× speedup for atlas generation. The core algorithm's correctness, time complexity, and efficiency-accuracy trade-offs are formally guaranteed using modern distance geometry, geometric constraint systems and combinatorial rigidity. The methodology further links the shape of the input assembling units to a type of intuitive and queryable bar-code of the output atlas, which in turn determine stable assembled structures and kinetics. This succinct input-output relationship facilitates reverse analysis and control toward design. A novel feature that is crucial to both the high sampling efficiency and decoupling of roadmap generation from sampling is a recently developed theory of convex Cayley (distance-based) custom parametrizations specific to assembly, as opposed to folding. Representing microstates with macrostate-specific Cayley parameters, to generate microstate samples, avoids gradient-descent search used by all prevailing methods. Further, these parametrizations convexify conformational regions or macrostates. This ratchets up sampling efficiency, significantly reducing number of repeated and discarded samples. These features of the new stand-alone methodology can also be used to complement the strengths of prevailing methodologies including Molecular Dynamics, Monte Carlo, and Fast Fourier Transform based methods.

The Role of Conformational Dynamics and Allostery in Modulating Protein Evolution

Advances in sequencing techniques and statistical methods have made it possible not only to predict sequences of ancestral proteins but also to identify thousands of mutations in the human exome, some of which are disease associated. These developments have motivated numerous theories and raised many questions regarding the fundamental principles behind protein evolution, which have been traditionally investigated horizontally using the tip of the phylogenetic tree through comparative studies of extant proteins within a family. In this article, we review a vertical comparison of the modern and resurrected ancestral proteins. We focus mainly on the dynamical properties responsible for a protein's ability to adapt new functions in response to environmental changes. Using the Dynamic Flexibility Index and the Dynamic Coupling Index to quantify the relative flexibility and dynamic coupling at a site-specific, single-amino-acid level, we provide evidence that the migration of hinges, which are often functionally critical rigid sites, is a mechanism through which proteins can rapidly evolve. Additionally, we show that disease-associated mutations in proteins often result in flexibility changes even at positions distal from mutational sites, particularly in the modulation of active site dynamics.

Combined strategies in structure-based virtual screening

The identification and optimization of lead compounds are inalienable components in drug design and discovery pipelines. As a powerful computational approach for the identification of hits with novel structural scaffolds, structure-based virtual screening (SBVS) has exhibited a remarkably increasing influence in the early stages of drug discovery. During the past decade, a variety of techniques and algorithms have been proposed and tested with different purposes in the scope of SBVS. Although SBVS has been a common and proven technology, it still shows some challenges and problems that are needed to be addressed, where the negative influence regardless of protein flexibility and the inaccurate prediction of binding affinity are the two major challenges. Here, focusing on these difficulties, we summarize a series of combined strategies or workflows developed by our group and others. Furthermore, several representative successful applications from recent publications are also discussed to demonstrate the effectiveness of the combined SBVS strategies in drug discovery campaigns.

CB-Dock: a web server for cavity detection-guided protein-ligand blind docking

As the number of elucidated protein structures is rapidly increasing, the growing data call for methods to efficiently exploit the structural information for biological and pharmaceutical purposes. Given the three-dimensional (3D) structure of a protein and a ligand, predicting their binding sites and affinity are a key task for computer-aided drug discovery. To address this task, a variety of docking tools have been developed. Most of them focus on docking in the preset binding sites given by users. To automatically predict binding modes without information about binding sites, we developed a user-friendly blind docking web server, named CB-Dock, which predicts binding sites of a given protein and calculates the centers and sizes with a novel curvature-based cavity detection approach, and performs docking with a popular docking program, Autodock Vina. This method was carefully optimized and achieved ~70% success rate for the top-ranking poses whose root mean square deviation (RMSD) were within 2 Å from the X-ray pose, which outperformed the state-of-the-art blind docking tools in our benchmark tests. CB-Dock offers an interactive 3D visualization of results, and is freely available at http://cao.labshare.cn/cb-dock/.

Holo-like and Druggable Protein Conformations from Enhanced Sampling of Binding Pocket Volume and Shape

Understanding molecular recognition of small molecules by proteins in atomistic detail is key for drug design. Molecular docking is a widely used computational method to mimic ligand-protein association in silico. However, predicting conformational changes occurring in proteins upon ligand binding is still a major challenge. Ensemble docking approaches address this issue by considering a set of different conformations of the protein obtained either experimentally or from computer simulations, e.g., molecular dynamics. However, holo structures prone to host (the correct) ligands are generally poorly sampled by standard molecular dynamics simulations of the apo protein. In order to address this limitation, we introduce a computational approach based on metadynamics simulations called ensemble docking with enhanced sampling of pocket shape (EDES) that allows holo-like conformations of proteins to be generated by exploiting only their apo structures. This is achieved by defining a set of collective variables that effectively sample different shapes of the binding site, ultimately mimicking the steric effect due to the ligand. We assessed the method on three challenging proteins undergoing different extents of conformational changes upon ligand binding. In all cases our protocol generates a significant fraction of structures featuring a low RMSD from the experimental holo geometry. Moreover, ensemble docking calculations using those conformations yielded in all cases native-like poses among the top-ranked ones.

Ligand discrimination between active and inactive activation loop conformations of Aurora-A kinase is unmodified by phosphorylation

Activation loop phosphorylation changes the position of equilibrium between DFG-in-like and DFG-out-like conformations but not the conformational preference of inhibitors.

Structure-based drug design is commonly used to guide the development of potent and specific enzyme inhibitors. Many enzymes – such as protein kinases – adopt multiple conformations, and conformational interconversion is expected to impact on the design of small molecule inhibitors. We measured the dynamic equilibrium between DFG-in-like active and DFG-out-like inactive conformations of the activation loop of unphosphorylated Aurora-A alone, in the presence of the activator TPX2, and in the presence of kinase inhibitors. The unphosphorylated kinase had a shorter residence time of the activation loop in the active conformation and a shift in the position of equilibrium towards the inactive conformation compared with phosphorylated kinase for all conditions measured. Ligand binding was associated with a change in the position of conformational equilibrium which was specific to each ligand and independent of the kinase phosphorylation state. As a consequence of this, the ability of a ligand to discriminate between active and inactive activation loop conformations was also independent of phosphorylation. Importantly, we discovered that the presence of multiple enzyme conformations can lead to a plateau in the overall ligand K_d, despite increasing affinity for the chosen target conformation, and modelled the conformational discrimination necessary for a conformation-promoting ligand.

Unraveling the Motions behind Enterovirus 71 Uncoating

Enterovirus 71 can be a severe pathogen in small children and immunocompromised adults. Virus uncoating is a critical step in the infection of the host cell; however, the mechanisms that control this process remain poorly understood. We applied normal mode analysis and perturbation response scanning to several complexes of the virus capsid and present a coarse-graining approach to analyze the full capsid. We show that our method offers an alternative to expressing the system as a set of rigid blocks and accounts for the interconnection between nodes within each subunit and protein interfaces across the capsid. In our coarse-grained approach, the modes associated with capsid expansion are captured in the first three nondegenerate modes and correspond to the changes observed in structural studies of the virus. We show that the resolution of the analysis may be modified without losing information on the global motions leading to uncoating. Perturbation response scanning revealed that a protomer cannot serve as a functional unit to explain deformations of the capsid. Instead, we define a pentamer as the minimum functional unit to investigate changes within the capsid. From the modal analysis and perturbation response scanning, we locate a hotspot region surrounding the fivefold axis. The range of the effect of these single, hotspot residues extend to 140 Å. The perturbation of internal capsid residues in this region displayed greatest propensity to capsid expansion, thus indicating the significant role that the RNA genome may play in triggering uncoating.

Ancient thioredoxins evolved to modern-day stability-function requirement by altering native state ensemble

Thioredoxins (THRXs)-small globular proteins that reduce other proteins-are ubiquitous in all forms of life, from Archaea to mammals. Although ancestral thioredoxins share sequential and structural similarity with the modern-day (extant) homologues, they exhibit significantly different functional activity and stability. We investigate this puzzle by comparative studies of their (ancient and modern-day THRXs') native state ensemble, as quantified by the dynamic flexibility index (DFI), a metric for the relative resilience of an amino acid to perturbations in the rest of the protein. Clustering proteins using DFI profiles strongly resemble an alternative classification scheme based on their activity and stability. The DFI profiles of the extant proteins are substantially different around the α3, α4 helices and catalytic regions. Likewise, allosteric coupling of the active site with the rest of the protein is different between ancient and extant THRXs, possibly explaining the decreased catalytic activity at low pH with evolution. At a global level, we note that the population of low-flexibility (called hinges) and high-flexibility sites increases with evolution. The heterogeneity (quantified by the variance) in DFI distribution increases with the decrease in the melting temperature typically associated with the evolution of ancient proteins to their modern-day counterparts.This article is part of a discussion meeting issue 'Allostery and molecular machines'.