A workflow for exploring ligand dissociation from a macromolecule: Efficient random acceleration molecular dynamics simulation and interaction fingerprint analysis of ligand trajectories

Multiresolution molecular dynamics simulations reveal the interplay between conformational variability and functional interactions in membrane-bound cytochrome P450 2B4

Cytochrome P450 2B4 (CYP 2B4) is one of the best-characterized CYPs and serves as a key model system for understanding the mechanisms of microsomal class II CYPs, which metabolize most known drugs. The highly flexible nature of CYP 2B4 is apparent from crystal structures that show the active site with either a wide open or a closed heme binding cavity. Here, we investigated the conformational ensemble of the full-length CYP 2B4 in a phospholipid bilayer, using multiresolution molecular dynamics (MD) simulations. Coarse-grained MD simulations revealed two predominant orientations of CYP 2B4's globular domain with respect to the bilayer. Their refinement by atomistic resolution MD showed adaptation of the enzyme's interaction with the lipid bilayer, leading to open configurations that facilitate ligand access to the heme binding cavity. CAVER analysis of enzyme tunnels, AquaDuct analysis of water routes, and Random Acceleration Molecular Dynamics simulations of ligand dissociation support the conformation-dependent passage of molecules between the active site and the protein surroundings. Furthermore, simulation of the re-entry of the inhibitor bifonazole into the open conformation of CYP 2B4 resulted in binding at a transient hydrophobic pocket within the active site cavity that may play a role in substrate binding or allosteric regulation. Together, these results show how the open conformation of CYP 2B4 facilitates the binding of substrates from and release of products to the membrane, whereas the closed conformation prolongs the residence time of substrates or inhibitors and selectively allows the passage of smaller reactants via the solvent and water channels.

Mining unique cysteine synthetases and computational study on thoroughly eliminating feedback inhibition through tunnel engineering

L-cysteine is an essential component in pharmaceutical and agricultural industries, and synthetic biology has made strides in developing new metabolic pathways for its production, particularly in archaea with unique O-phosphoserine sulfhydrylases (OPSS) as key enzymes. In this study, we employed database mining to identify a highly catalytic activity OPSS from Acetobacterium sp. (AsOPSS). However, it was observed that the enzymatic activity of AsOPSS suffered significant feedback inhibition from the product L-cysteine, exhibiting an IC50 value of merely 1.2 mM. A semi-rational design combined with tunnel analysis strategy was conducted to engineer AsOPSS. The best variant, AsOPSSA218R was achieved, totally eliminating product inhibition without sacrificing catalytic efficiency. Molecular docking and molecular dynamic simulations indicated that the binding conformation of AsOPSSA218R with L-cys was altered, leading to a reduced affinity between L-cysteine and the active pocket. Tunnel analysis revealed that the AsOPSSA218R variant reshaped the landscape of the tunnel, resulting in the construction of a new tunnel. Furthermore, random acceleration molecular dynamics simulation and umbrella sampling simulation demonstrated that the novel tunnel improved the suitability for product release and effectively separated the interference between the product release and substrate binding processes. Finally, more than 45 mM of L-cysteine was produced in vitro within 2 h using the AsOPSSA218R variant. Our findings emphasize the potential for relieving feedback inhibition by artificially generating new product release channels, while also laying an enzymatic foundation for efficient L-cysteine production.

Computational screening of the effects of mutations on protein-protein off-rates and dissociation mechanisms by τRAMD

The dissociation rate, or its reciprocal, the residence time (τ), is a crucial parameter for understanding the duration and biological impact of biomolecular interactions. Accurate prediction of τ is essential for understanding protein-protein interactions (PPIs) and identifying potential drug targets or modulators for tackling diseases. Conventional molecular dynamics simulation techniques are inherently constrained by their limited timescales, making it challenging to estimate residence times, which typically range from minutes to hours. Building upon its successful application in protein-small molecule systems, τ-Random Acceleration Molecular Dynamics (τRAMD) is here investigated for estimating dissociation rates of protein-protein complexes. τRAMD enables the observation of unbinding events on the nanosecond timescale, facilitating rapid and efficient computation of relative residence times. We tested this methodology for three protein-protein complexes and their extensive mutant datasets, achieving good agreement between computed and experimental data. By combining τRAMD with MD-IFP (Interaction Fingerprint) analysis, dissociation mechanisms were characterized and their sensitivity to mutations investigated, enabling the identification of molecular hotspots for selective modulation of dissociation kinetics. In conclusion, our findings underscore the versatility of τRAMD as a simple and computationally efficient approach for computing relative protein-protein dissociation rates and investigating dissociation mechanisms, thereby aiding the design of PPI modulators.

Elucidating the Selective Mechanism of Drugs Targeting Cyclin-Dependent Kinases with Integrated MetaD-US Simulation

Cyclin-dependent kinases (CDKs), including CDK12 and CDK13, play crucial roles in regulating the cell cycle and RNA polymerase II activity, making them vital targets for cancer therapies. SR4835 is a selective inhibitor of CDK12/13, showing significant potential for treating triple-negative breast cancer. To elucidate the selective mechanism of SR4835 among three CDKs (CDK13/12/9), we developed an innovative enhanced sampling method, integrated well-tempered metadynamics-umbrella sampling (IMUS). IMUS synergistically combines the comprehensive pathway exploration capability of well-tempered metadynamics (WT-MetaD) with the precise free energy calculation capability of umbrella sampling, enabling the efficient and accurate characterization of drug-target interactions. The accurate calculation of binding free energy and the detailed analysis of the kinetic mechanism of the drug-target interaction using IMUS successfully elucidate the drug selectivity mechanism targeting the three CDKs, showing that the selectivity is primarily arising from differences in the stability of H-bonds within the Hinge region of the kinases and the interaction patterns during the protein-ligand recognition process. These findings also underscore the utility of IMUS in efficiently and accurately capturing drug-target interaction processes with clear mechanisms.

Methods for Classical-Mechanical Molecular Simulation in Chemistry: Achievements, Limitations, Perspectives

More than a half century ago it became feasible to simulate, using classical-mechanical equations of motion, the dynamics of molecular systems on a computer. Since then classical-physical molecular simulation has become an integral part of chemical research. It is widely applied in a variety of branches of chemistry and has significantly contributed to the development of chemical knowledge. It offers understanding and interpretation of experimental results, semiquantitative predictions for measurable and nonmeasurable properties of substances, and allows the calculation of properties of molecular systems under conditions that are experimentally inaccessible. Yet, molecular simulation is built on a number of assumptions, approximations, and simplifications which limit its range of applicability and its accuracy. These concern the potential-energy function used, adequate sampling of the vast statistical-mechanical configurational space of a molecular system and the methods used to compute particular properties of chemical systems from statistical-mechanical ensembles. During the past half century various methodological ideas to improve the efficiency and accuracy of classical-physical molecular simulation have been proposed, investigated, evaluated, implemented in general simulation software or were abandoned. The latter because of fundamental flaws or, while being physically sound, computational inefficiency. Some of these methodological ideas are briefly reviewed and the most effective methods are highlighted. Limitations of classical-physical simulation are discussed and perspectives are sketched.

6-Phosphogluconate dehydrogenase 2 bridges the OPP and shikimate pathways to enhance aromatic amino acid production in plants

The oxidative pentose phosphate (OPP) pathway provides metabolic intermediates for the shikimate pathway and directs carbon flow to the biosynthesis of aromatic amino acids (AAAs), which serve as basic protein building blocks and precursors of numerous metabolites essential for plant growth. However, genetic evidence linking the two pathways is largely unclear. In this study, we identified 6-phosphogluconate dehydrogenase 2 (PGD2), the rate-limiting enzyme of the cytosolic OPP pathway, through suppressor screening of arogenate dehydrogenase 2 (adh2) in Arabidopsis. Our data indicated that a single amino acid substitution at position 63 (glutamic acid to lysine) of PGD2 enhanced its enzyme activity by facilitating the dissociation of products from the active site of PGD2, thus increasing the accumulation of AAAs and partially restoring the defective phenotype of adh2. Phylogenetic analysis indicated that the point mutation occurred in a well-conserved amino acid residue. Plants with different amino acids at this conserved site of PGDs confer diverse catalytic activities, thus exhibiting distinct AAAs producing capability. These findings uncover the genetic link between the OPP pathway and AAAs biosynthesis through PGD2. The gain-of-function point mutation of PGD2 identified here could be considered as a potential engineering target to alter the metabolic flux for the production of AAAs and downstream compounds.

Structure-guided discovery of novel AflG inhibitors for aflatoxin contamination control in aspergillus flavus

Aflatoxins (AFs) are highly carcinogenic metabolites produced by Aspergillus species that can contaminate critical food staples, leading to significant health and economic risks. The cytochrome P450 monooxygenase AflG catalyzes an early step in AF biosynthesis, resulting in the conversion of averantin (AVN) to 5'-hydroxy-averantin. However, the molecular mechanism underlying the AflG-AVN interaction remains unclear. Here, we sought to understand the structural features of AflG in complex with AVN to enable the identification of inhibitors targeting the AflG binding pocket. To achieve this goal, we employed a comprehensive approach combining computational and experimental methods. Structural modeling and microsecond-scale molecular dynamics (MD) simulations yielded new insights into AflG architecture and unveiled unique ligand binding conformations of the AflG-AVN complex. High-throughput virtual screening of more than 1.3 million compounds pinpointed specific subsets with favorable predicted docking scores. The resulting compounds were ranked based on binding free energy calculations and evaluated with MD simulations and in vitro experiments with Aspergillus flavus. Our results revealed two compounds significantly inhibited AF biosynthesis. Comprehensive structural analysis elucidated the binding sites of competitive inhibitors and demonstrated their regulation of AflG dynamics. This structure-guided pipeline successfully enabled the identification of novel AflG inhibitors and provided novel molecular insights that will guide future efforts to develop effective therapeutics that prevent AF contamination.

Deciphering the Selectivity of CBL-B Inhibitors Using All-Atom Molecular Dynamics and Machine Learning

We employ a combination of accelerated molecular dynamics and machine learning to unravel how the dynamic characteristics of CBL-B and C-CBL confer their binding affinity and selectivity for ligands from subtle structural disparities within their binding pockets and dissociation pathways. Our predictive model of dissociation rate constants (k off) demonstrates a moderate correlation between predicted k off and experimental IC50 values, which is consistent with experimental k off and τ-random accelerated molecular dynamics (τRAMD) results. By employing a linear regression of dissociation trajectories, we identified key amino acids in binding pockets and along the dissociation paths responsible for activity and selectivity. These amino acids are statistically significant in achieving activity and selectivity and contribute to the primary structural discrepancies between CBL-B and C-CBL. Moreover, the binding free energies calculated from molecular mechanics with generalized Born and surface area solvation (MM/GBSA) highlight the ΔG difference between CBL-B and C-CBL. The k off prediction, together with the key amino acids, provides important guides for designing drugs with high selectivity.

Prebound State Discovered in the Unbinding Pathway of Fluorinated Variants of the Trypsin-BPTI Complex Using Random Acceleration Molecular Dynamics Simulations

The serine protease trypsin forms a tightly bound inhibitor complex with the bovine pancreatic trypsin inhibitor (BPTI). The complex is stabilized by the P1 residue Lys15, which interacts with negatively charged amino acids at the bottom of the S1 pocket. Truncating the P1 residue of wildtype BPTI to α-aminobutyric acid (Abu) leaves a complex with moderate inhibitor strength, which is held in place by additional hydrogen bonds at the protein-protein interface. Fluorination of the Abu residue partially restores the inhibitor strength. The mechanism with which fluorination can restore the inhibitor strength is unknown, and accurate computational investigation requires knowledge of the binding and unbinding pathways. The preferred unbinding pathway is likely to be complex, as encounter states have been described before, and unrestrained umbrella sampling simulations of these complexes suggest additional energetic minima. Here, we use random acceleration molecular dynamics to find a new metastable state in the unbinding pathway of Abu-BPTI variants and wildtype BPTI from trypsin, which we call the prebound state. The prebound state and the fully bound state differ by a substantial shift in the position, a slight shift in the orientation of the BPTI variants, and changes in the interaction pattern. Particularly important is the breaking of three hydrogen bonds around Arg17. Fluorination of the P1 residue lowers the energy barrier of the transition between the fully bound state and prebound state and also lowers the energy minimum of the prebound state. While the effect of fluorination is in general difficult to quantify, here, it is in part caused by favorable stabilization of a hydrogen bond between Gln194 and Cys14. The interaction pattern of the prebound state offers insights into the inhibitory mechanism of BPTI and might add valuable information for the design of serine protease inhibitors.

Learning Protein-Ligand Unbinding Pathways via Single-Parameter Community Detection

Understanding the dynamics of biomolecular complexes, e.g., of protein-ligand (un)binding, requires the comprehension of paths such systems take between metastable states. In MD simulations, paths are usually not observable per se, but they need to be inferred from simulation trajectories. Here, we present a novel approach to cluster trajectories based on a community detection algorithm that necessitates only the definition of a single parameter. The unbinding of the streptavidin-biotin complex is used as a benchmark system and the A2a adenosine receptor in complex with the inhibitor ZM241385 as an elaborate application. We demonstrate how such clusters of trajectories correspond to pathways and how the approach helps in the identification of reaction coordinates for a considered (un)binding process.

Understanding the Differences of Danusertib's Residence Time in Aurora Kinases A/B: Dissociation Paths and Key Residues Identified using Conventional and Enhanced Molecular Dynamics Simulations

The accurate experimental estimation of protein-ligand systems' residence time (τ) has become very relevant in drug design projects due to its importance in the last stages of refinement of the drug's pharmacodynamics and pharmacokinetics. It is now well-known that it is not sufficient to estimate the affinity of a protein-drug complex in the thermodynamic equilibrium process in in vitro experiments (closed systems), where the concentrations of the drug and protein remain constant. On the contrary, it is mandatory to consider the conformational dynamics of the system in terms of the binding and unbinding processes between protein and drugs in in vivo experiments (open systems), where their concentrations are in constant flux. This last model has been proven to dictate much of several drugs' pharmacological activities in vivo. At the atomistic level, molecular dynamics simulations can explain why some drugs are more effective than others or unveil the molecular aspects that make some drugs work better in one molecular target. Here, the protein kinases Aurora A/B, complexed with its inhibitor Danusertib, were studied using conventional and enhanced molecular dynamics (MD) simulations to estimate the dissociation paths and, therefore, the computational τ values and their comparison with experimental ones. Using classical molecular dynamics (cMD), three differential residues within the Aurora A/B active site, which seems to play an essential role in the observed experimental Danusertib's residence time against these kinases, were characterized. Then, using WT-MetaD, the relative Danusertib's residence times against Aurora A/B kinases were measured in a nanosecond time scale and were compared to those τ values observed experimentally. In addition, the potential dissociation paths of Danusertib in Aurora A and B were characterized, and differences that might be explained by the differential residues in the enzyme's active sites were found. In perspective, it is expected that this computational protocol can be applied to other protein-ligand complexes to understand, at the molecular level, the differences in residence times and amino acids that may contribute to it.

Dynamical Reweighting for Biased Rare Event Simulations

Dynamical reweighting techniques aim to recover the correct molecular dynamics from a simulation at a modified potential energy surface. They are important for unbiasing enhanced sampling simulations of molecular rare events. Here, we review the theoretical frameworks of dynamical reweighting for modified potentials. Based on an overview of kinetic models with increasing level of detail, we discuss techniques to reweight two-state dynamics, multistate dynamics, and path integrals. We explore the natural link to transition path sampling and how the effect of nonequilibrium forces can be reweighted. We end by providing an outlook on how dynamical reweighting integrates with techniques for optimizing collective variables and with modern potential energy surfaces.

Characterization of the Bottlenecks and Pathways for Inhibitor Dissociation from [NiFe] Hydrogenase

[NiFe] hydrogenases can act as efficient catalysts for hydrogen oxidation and biofuel production. However, some [NiFe] hydrogenases are inhibited by gas molecules present in the environment, such as O2 and CO. One strategy to engineer [NiFe] hydrogenases and achieve O2- and CO-tolerant enzymes is by introducing point mutations to block the access of inhibitors to the catalytic site. In this work, we characterized the unbinding pathways of CO in the complex with the wild-type and 10 different mutants of [NiFe] hydrogenase from Desulfovibrio fructosovorans using τ-random accelerated molecular dynamics (τRAMD) to enhance the sampling of unbinding events. The ranking provided by the relative residence times computed with τRAMD is in agreement with experiments. Extensive data analysis of the simulations revealed that from the two bottlenecks proposed in previous studies for the transit of gas molecules (residues 74 and 122 and residues 74 and 476), only one of them (residues 74 and 122) effectively modulates diffusion and residence times for CO. We also computed pathway probabilities for the unbinding of CO, O2, and H2 from the wild-type [NiFe] hydrogenase, and we observed that while the most probable pathways are the same, the secondary pathways are different. We propose that introducing mutations to block the most probable paths, in combination with mutations to open the main secondary path used by H2, can be a feasible strategy to achieve CO and O2 resistance in the [NiFe] hydrogenase from Desulfovibrio fructosovorans.

Role of UDP-N-acetylmuramic acid in the regulation of MurA activity revealed by molecular dynamics simulations

The peptidoglycan biosynthesis pathway plays a vital role in bacterial cells, and facilitates peptidoglycan layer formation, a fundamental structural component of the bacterial cell wall. The enzymes in this pathway are candidates for antibiotic development, as most do not have mammalian homologues. The UDP-N-acetylglucosamine (UNAG) enolpyruvyl transferase enzyme (MurA) in the peptidoglycan pathway cytoplasmic step is responsible for the phosphoenolpyruvate (PEP)-UNAG catalytic reaction, forming UNAG enolpyruvate and inorganic phosphate. Reportedly, UDP-N-acetylmuramic acid (UNAM) binds tightly to MurA forming a dormant UNAM-PEP-MurA complex and acting as a MurA feedback inhibitor. MurA inhibitors are complex, owing to competitive binding interactions with PEP, UNAM, and UNAG at the MurA active site. We used computational methods to explore UNAM and UNAG binding. UNAM showed stronger hydrogen-bond interactions with the Arg120 and Arg91 residues, which help to stabilize the closed conformation of MurA, than UNAG. Binding free energy calculations using end-point computational methods showed that UNAM has a higher binding affinity than UNAG, when PEP is attached to Cys115. The unbinding process, simulated using τ-random acceleration molecular dynamics, showed that UNAM has a longer relative residence time than UNAG, which is related to several complex dissociation pathways, each with multiple intermediate metastable states. This prevents the loop from opening and exposing the Arg120 residue to accommodate UNAG and potential new ligands. Moreover, we demonstrate the importance of Cys115-linked PEP in closed-state loop stabilization. We provide a basis for evaluating novel UNAM analogues as potential MurA inhibitors. PUBLIC SIGNIFICANCE: MurA is a critical enzyme involved in bacterial cell wall biosynthesis and is involved in antibiotic resistance development. UNAM can remain in the target protein's active site for an extended time compared to its natural substrate, UNAG. The prolonged interaction of this highly stable complex known as the 'dormant complex' comprises UNAM-PEP-MurA and offers insights into antibiotic development, providing potential options against drug-resistant bacteria and advancing our understanding of microbial biology.

From Molecular Dynamics to Taste Sensory Perception: A Comprehensive Study on the Interaction of Umami Peptides with the T1R1/T1R3-VFT Receptor

The research on the umami receptor-ligand interaction is crucial for understanding umami perception. This study integrated molecular simulations, sensory evaluation, and biosensor technology to analyze the interaction between umami peptides and the umami receptor T1R1/T1R3-VFT. Molecular dynamics simulations were used to investigate the dissociation process of seven umami peptides with the umami receptor T1R1/T1R3-VFT, and by calculating the potential mean force curve using the Jarzynski equation, it was found that the binding free energy of umami peptide is between -58.80 and -12.17 kcal/mol, which had a strong correlation with the umami intensity obtained by time intensity sensory evaluation. Through correlation analysis, the dissociation rate constants (0.0126-0.394 1/s) of umami peptides were found to have a great impact on umami perception. The faster the dissociation rate of umami peptides from receptors, the stronger the perceived intensity of the umami taste. This research aims to elucidate the relationship between the umami peptide-receptor interaction and umami perception, providing theoretical support for the exploration of umami perception mechanisms.

Systematic analysis, aggregation and visualisation of interaction fingerprints for molecular dynamics simulation data

Computational methods such as molecular docking or molecular dynamics (MD) simulations have been developed to simulate and explore the interactions between biomolecules. However, the interactions obtained using these methods are difficult to analyse and evaluate. Interaction fingerprints (IFPs) have been proposed to derive interactions from static 3D coordinates and transform them into 1D bit vectors. More recently, the concept has been applied to derive IFPs from MD simulations, which adds a layer of complexity by adding the temporal motion and dynamics of a system. As a result, many IFPs are obtained from one MD simulation, resulting in a large number of individual IFPs that are difficult to analyse compared to IFPs derived from static 3D structures. Scientific contribution: We introduce a new method to systematically aggregate IFPs derived from MD simulation data. In addition, we propose visualisations to effectively analyse and compare IFPs derived from MD simulation data to account for the temporal evolution of interactions and to compare IFPs across different MD simulations. This has been implemented as a freely available Python library and can therefore be easily adopted by other researchers and to different MD simulation datasets.

Computational Workflow for Refining AlphaFold Models in Drug Design Using Kinetic and Thermodynamic Binding Calculations: A Case Study for the Unresolved Inactive Human Adenosine A3 Receptor

A structure-based drug design pipeline that considers both thermodynamic and kinetic binding data of ligands against a receptor will enable the computational design of improved drug molecules. For unresolved GPCR-ligand complexes, a workflow that can apply both thermodynamic and kinetic binding data in combination with alpha-fold (AF)-derived or other homology models and experimentally resolved binding modes of relevant ligands in GPCR-homologs needs to be tested. Here, as test case, we studied a congeneric set of ligands that bind to a structurally unresolved G protein-coupled receptor (GPCR), the inactive human adenosine A3 receptor (hA3R). We tested three available homology models from which two have been generated from experimental structures of hA1R or hA2AR and one model was a multistate alphafold 2 (AF2)-derived model. We applied alchemical calculations with thermodynamic integration coupled with molecular dynamics (TI/MD) simulations to calculate the experimental relative binding free energies and residence time (τ)-random accelerated MD (τ-RAMD) simulations to calculate the relative residence times (RTs) for antagonists. While the TI/MD calculations produced, for the three homology models, good Pearson correlation coefficients, correspondingly, r = 0.74, 0.62, and 0.67 and mean unsigned error (mue) values of 0.94, 1.31, and 0.81 kcal mol-1, the τ-RAMD method showed r = 0.92 and 0.52 for the first two models but failed to produce accurate results for the multistate AF2-derived model. With subsequent optimization of the AF2-derived model by reorientation of the side chain of R1735.34 located in the extracellular loop 2 (EL2) that blocked ligand's unbinding, the computational model showed r = 0.84 for kinetic data and improved performance for thermodynamic data (r = 0.81, mue = 0.56 kcal mol-1). Overall, after refining the multistate AF2 model with physics-based tools, we were able to show a strong correlation between predicted and experimental ligand relative residence times and affinities, achieving a level of accuracy comparable to an experimental structure. The computational workflow used can be applied to other receptors, helping to rank candidate drugs in a congeneric series and enabling the prioritization of leads with stronger binding affinities and longer residence times.

Fingerprinting Interactions between Proteins and Ligands for Facilitating Machine Learning in Drug Discovery

Molecular recognition is fundamental in biology, underpinning intricate processes through specific protein-ligand interactions. This understanding is pivotal in drug discovery, yet traditional experimental methods face limitations in exploring the vast chemical space. Computational approaches, notably quantitative structure-activity/property relationship analysis, have gained prominence. Molecular fingerprints encode molecular structures and serve as property profiles, which are essential in drug discovery. While two-dimensional (2D) fingerprints are commonly used, three-dimensional (3D) structural interaction fingerprints offer enhanced structural features specific to target proteins. Machine learning models trained on interaction fingerprints enable precise binding prediction. Recent focus has shifted to structure-based predictive modeling, with machine-learning scoring functions excelling due to feature engineering guided by key interactions. Notably, 3D interaction fingerprints are gaining ground due to their robustness. Various structural interaction fingerprints have been developed and used in drug discovery, each with unique capabilities. This review recapitulates the developed structural interaction fingerprints and provides two case studies to illustrate the power of interaction fingerprint-driven machine learning. The first elucidates structure-activity relationships in β2 adrenoceptor ligands, demonstrating the ability to differentiate agonists and antagonists. The second employs a retrosynthesis-based pre-trained molecular representation to predict protein-ligand dissociation rates, offering insights into binding kinetics. Despite remarkable progress, challenges persist in interpreting complex machine learning models built on 3D fingerprints, emphasizing the need for strategies to make predictions interpretable. Binding site plasticity and induced fit effects pose additional complexities. Interaction fingerprints are promising but require continued research to harness their full potential.

How Robust Is the Ligand Binding Transition State?

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 binding transition state is possible with the help of enhanced sampling molecular dynamics methods. Previously, we generated unbinding paths for an inhibitor of soluble epoxide hydrolase (sEH) with a residence time of 11 min. Here, we computationally modeled unbinding events with the weighted ensemble method REVO (resampling of ensembles by variation optimization) for five additional inhibitors of sEH with residence times ranging from 14.25 to 31.75 min, with average prediction accuracy within an order of magnitude. The unbinding ensembles are analyzed in detail, focusing on features of the ligand binding transition state ensembles (TSEs). We find that ligands with similar bound poses can show significant differences in their ligand binding TSEs, in terms of their spatial distribution and protein-ligand interactions. However, we also find similarities across the TSEs when examining more general features such as ligand degrees of freedom. Together these findings show significant challenges for rational, kinetics-based drug design.

Structural Insights into Dopamine Receptor-Ligand Interactions: From Agonists to Antagonists

This study explores the intricacies of dopamine receptor-ligand interactions, focusing on the D1R and D5R subtypes. Using molecular modeling techniques, we investigate the binding of the pan-agonist rotigotine, revealing a universal binding mode at the orthosteric binding pocket (OBP). Additionally, we analyze the stability of antagonist-receptor complexes with SKF83566 and SCH23390. By examining the impact of specific mutations on ligand-receptor interactions through computational simulations and thermostability assays, we gain insights into binding stability. Our research also delves into the structural and energetic aspects of antagonist binding to D1R and D5R in their inactive states. These findings enhance our understanding of dopamine receptor pharmacology and hold promise for drug development in central nervous system disorders, opening doors to future research and innovation in this field.

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.

MD-Ligand-Receptor: A High-Performance Computing Tool for Characterizing Ligand-Receptor Binding Interactions in Molecular Dynamics Trajectories

Molecular dynamics simulation is a widely employed computational technique for studying the dynamic behavior of molecular systems over time. By simulating macromolecular biological systems consisting of a drug, a receptor and a solvated environment with thousands of water molecules, MD allows for realistic ligand-receptor binding interactions (lrbi) to be studied. In this study, we present MD-ligand-receptor (MDLR), a state-of-the-art software designed to explore the intricate interactions between ligands and receptors over time using molecular dynamics trajectories. Unlike traditional static analysis tools, MDLR goes beyond simply taking a snapshot of ligand-receptor binding interactions (lrbi), uncovering long-lasting molecular interactions and predicting the time-dependent inhibitory activity of specific drugs. With MDLR, researchers can gain insights into the dynamic behavior of complex ligand-receptor systems. Our pipeline is optimized for high-performance computing, capable of efficiently processing vast molecular dynamics trajectories on multicore Linux servers or even multinode HPC clusters. In the latter case, MDLR allows the user to analyze large trajectories in a very short time. To facilitate the exploration and visualization of lrbi, we provide an intuitive Python notebook (Jupyter), which allows users to examine and interpret the results through various graphical representations.

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.

Structure-Kinetic Relationship for Drug Design Revealed by a PLS Model with Retrosynthesis-Based Pre-Trained Molecular Representation and Molecular Dynamics Simulation

Drug design based on kinetic properties is growing in application. Here, we applied retrosynthesis-based pre-trained molecular representation (RPM) in machine learning (ML) to train 501 inhibitors of 55 proteins and successfully predicted the dissociation rate constant (k_off) values of 38 inhibitors from an independent dataset for the N-terminal domain of heat shock protein 90α (N-HSP90). Our RPM molecular representation outperforms other pre-trained molecular representations such as GEM, MPG, and general molecular descriptors from RDKit. Furthermore, we optimized the accelerated molecular dynamics to calculate the relative retention time (RT) for the 128 inhibitors of N-HSP90 and obtained the protein-ligand interaction fingerprints (IFPs) on their dissociation pathways and their influencing weights on the k_off value. We observed a high correlation among the simulated, predicted, and experimental -log(k_off) values. Combining ML, molecular dynamics (MD) simulation, and IFPs derived from accelerated MD helps design a drug for specific kinetic properties and selectivity profiles to the target of interest. To further validate our k_off predictive ML model, we tested our model on two new N-HSP90 inhibitors, which have experimental k_off values and are not in our ML training dataset. The predicted k_off values are consistent with experimental data, and the mechanism of their kinetic properties can be explained by IFPs, which shed light on the nature of their selectivity against N-HSP90 protein. We believe that the ML model described here is transferable to predict k_off of other proteins and will enhance the kinetics-based drug design endeavor.

Predicting Protein-Ligand Binding and Unbinding Kinetics with Biased MD Simulations and Coarse-Graining of Dynamics: Current State and Challenges

The prediction of drug-target binding and unbinding kinetics that occur on time scales between milliseconds and several hours is a prime challenge for biased molecular dynamics simulation approaches. This Perspective gives a concise summary of the theory and the current state-of-the-art of such predictions via biased simulations, of insights into the molecular mechanisms defining binding and unbinding kinetics as well as of the extraordinary challenges predictions of ligand kinetics pose in comparison to binding free energy predictions.

Novel Inhibitory Role of Fenofibric Acid by Targeting Cryptic Site on the RBD of SARS-CoV-2

The emergence of the recent pandemic causing severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has created an alarming situation worldwide. It also prompted extensive research on drug repurposing to find a potential treatment for SARS-CoV-2 infection. An active metabolite of the hyperlipidemic drug fenofibrate (also called fenofibric acid or FA) was found to destabilize the receptor-binding domain (RBD) of the viral spike protein and therefore inhibit its binding to human angiotensin-converting enzyme 2 (hACE2) receptor. Despite being considered as a potential drug candidate for SARS-CoV-2, FA's inhibitory mechanism remains to be elucidated. We used molecular dynamics (MD) simulations to investigate the binding of FA to the RBD of the SARS-CoV-2 spike protein and revealed a potential cryptic FA binding site. Free energy calculations were performed for different FA-bound RBD complexes. The results suggest that the interaction of FA with the cryptic binding site of RBD alters the conformation of the binding loop of RBD and effectively reduces its binding affinity towards ACE2. Our study provides new insights for the design of SARS-CoV-2 inhibitors targeting cryptic sites on the RBD of SARS-CoV-2.

Cumulative Millisecond-Long Sampling for a Comprehensive Energetic Evaluation of Aqueous Ionic Liquid Effects on Amino Acid Interactions

The interactions of amino acid side-chains confer diverse energetic contributions and physical properties to a protein's stability and function. Various computational tools estimate the effect of changing a given amino acid on the protein's stability based on parametrized (free) energy functions. When parametrized for the prediction of protein stability in water, such energy functions can lead to suboptimal results for other solvents, such as ionic liquids (IL), aqueous ionic liquids (aIL), or salt solutions. However, to our knowledge, no comprehensive data are available describing the energetic effects of aIL on intramolecular protein interactions. Here, we present the most comprehensive set of potential of mean force (PMF) profiles of pairwise protein-residue interactions to date, covering 50 relevant interactions in water, the two biotechnologically relevant aIL [BMIM/Cl] and [BMIM/TfO], and [Na/Cl]. These results are based on a cumulated simulation time of >1 ms. aIL and salt ions can weaken, but also strengthen, specific residue interactions by more than 3 kcal mol^-1, depending on the residue pair, residue-residue configuration, participating ions, and concentration, necessitating considering such interactions specifically. These changes originate from a complex interplay of competitive or cooperative noncovalent ion-residue interactions, changes in solvent structural dynamics, or unspecific charge screening effects and occur at the contact distance but also at larger, solvent-separated distances. This data provide explanations at the atomistic and energetic levels for complex IL effects on protein stability and should help improve the prediction accuracies of computational tools that estimate protein stability based on (free) energy functions.

Linolenic Acid-Derived Oxylipins Inhibit Aflatoxin Biosynthesis in Aspergillus flavus through Activation of Imizoquin Biosynthesis

Oxylipins play important signaling roles in aflatoxin (AF) biosynthesis in Aspergillus flavus. We previously showed that exogenous supply of autoxidated linolenic acid (AL) inhibited AF biosynthesis in A. flavus via oxylipins, but the molecular mechanism is still unknown. Here, we performed multiomics analyses of A. flavus grown in media with or without AL. Targeted metabolite analyses and quantitative reverse transcription (qRT)-polymerase chain reaction (PCR) showed that the imizoquin (IMQ) biosynthetic pathway was distinctly upregulated in the presence of AL. ¹³C-glucose labeling confirmed in parallel that the tricarboxylic acid cycle was also enhanced by AL, consistent with observed increases in mycelial growth. Moreover, we integrated thermal proteome profiling and molecular dynamics simulations to identify a potential receptor of AL; AL was found to interact with a transporter (ImqJ) located in the IMQ gene cluster, primarily through hydrophobic interactions. Further analyses of strains with an IMQ pathway transcription factor overexpressed or knocked out confirmed that this pathway was critical for AL-mediated inhibition of AF biosynthesis. Comparison of 22 assembled A. flavus and Aspergillus oryzae genomes showed that genes involved in the IMQ pathway were positively selected in A. oryzae. Taken together, the results of our study provide novel insights into oxylipin-mediated regulation of AF biosynthesis and suggest potential methods for preventing AF contamination of crops.

Free energy and kinetic rate calculation via non-equilibrium molecular simulation: application to biomolecules

Non-equilibrium molecular dynamics (NEMD) simulation has been recognized as a powerful tool for examining biomolecules and provides fruitful insights into not only non-equilibrium but also equilibrium processes. We review recent advances in NEMD simulation and relevant, fundamental results of non-equilibrium statistical mechanics. We first introduce Crooks fluctuation theorem and Jarzynski equality that relate free energy difference to work done on a physical system during a non-equilibrium process. The theorems are beneficial for the analysis of NEMD trajectories. We then describe rate theory, a framework to calculate molecular kinetics from a non-equilibrium process; this theoretical framework enables us to calculate a reaction time-mean-first passage time-from NEMD trajectories. We, in turn, present recent NEMD techniques that apply an external force to a system to enhance molecular dissociation and introduce their application to biomolecules. Lastly, we show the current status of an appropriate selection of reaction coordinates for NEMD simulation.

Bell-Evans model and steered molecular dynamics in uncovering the dissociation kinetics of ligands targeting G-protein-coupled receptors

Recently, academic and industrial scientific communities involved in kinetics-based drug development have become immensely interested in predicting the drug target residence time. Screening drug candidates in terms of their computationally predicted residence times, which is a measure of drug efficacy in vivo, and simultaneously assessing computational binding affinities are becoming inevitable. Non-equilibrium molecular simulation approaches are proven to be useful in this purpose. Here, we have implemented an optimized approach of combining the data derived from steered molecular dynamics simulations and the Bell-Evans model to predict the absolute residence times of the antagonist ZMA241385 and agonist NECA that target the A2A adenosine receptor of the G-protein-coupled receptor (GPCR) protein family. We have predicted the absolute ligand residence times on the timescale of seconds. However, our predictions were many folds shorter than those determined experimentally. Additionally, we calculated the thermodynamics of ligand binding in terms of ligand binding energies and the per-residue contribution of the receptor. Subsequently, binding pocket hotspot residues that would be important for further computational mutagenesis studies were identified. In the experiment, similar sets of residues were found to be in significant contact with both ligands under study. Our results build a strong foundation for further improvement of our approach by rationalizing the kinetics of ligand unbinding with the thermodynamics of ligand binding.

Protein Flexibility and Dissociation Pathway Differentiation Can Explain Onset of Resistance Mutations in Kinases

Understanding how mutations render a drug ineffective is a problem of immense relevance. Often the mechanism through which mutations cause drug resistance can be explained purely through thermodynamics. However, the more perplexing situation is when two proteins have the same drug binding affinities but different residence times. In this work, we demonstrate how all-atom molecular dynamics simulations using recent developments grounded in statistical mechanics can provide a detailed mechanistic rationale for such variances. We discover dissociation mechanisms for the anti-cancer drug Imatinib (Gleevec) against wild-type and the N368S mutant of Abl kinase. We show how this point mutation triggers far-reaching changes in the protein's flexibility and leads to a different, much faster, drug dissociation pathway. We believe that this work marks an efficient and scalable approach to obtain mechanistic insight into resistance mutations in biomolecular receptors that are hard to explain using a structural perspective.

Mechanistic Insight on General Protein-Binding Ability of ATP and the Impacts of Arginine Residues

Recent experiments suggested that adenosine triphosphate (ATP) can regulate liquid-liquid phase separation (LLPS) of various proteins and inhibit protein aggregations at its physiological concentration, which is highly correlated with the nonspecific interactions of ATP to a wide variety of proteins. However, the mechanism underlying the general binding capability of ATP largely remains unclear. In this work, we used molecular dynamics simulation to study the binding of ATPs to three proteins with distinct net charges: TDP-43 NTD (-7 e), TAF15-RRM (0 e), HWEL (+8 e). Negatively charged ATP exhibits a strong trend to accumulate around all of these proteins. While only a fraction of the accumulated ATPs directly binds to the limited regions of the protein surface, additional ATPs indirectly bind to proteins by aggregating into ATP clusters. Hence, the proportion of the directly bound ATPs in the clusters as well as their binding regions can be adjusted in response to different proteins, which makes ATP well adapted to a variety of proteins. Moreover, our results suggest that ATP tightly binds to Arg with high affinity, and Arg dominates the direct binding of ATP. Meanwhile, Arg also affects the self-association of accumulated ATPs. The size of the ATP cluster is effectively regulated by the distribution of Arg. Considering the ubiquity of Arg in proteins, our findings are helpful to understand the general binding capability of ATP.

Enhanced-Sampling Simulations for the Estimation of Ligand Binding Kinetics: Current Status and Perspective

The dissociation rate (k off) associated with ligand unbinding events from proteins is a parameter of fundamental importance in drug design. Here we review recent major advancements in molecular simulation methodologies for the prediction of k off. Next, we discuss the impact of the potential energy function models on the accuracy of calculated k off values. Finally, we provide a perspective from high-performance computing and machine learning which might help improve such predictions.

Determination of Multidirectional Pathways for Ligand Release from the Receptor: A New Approach Based on Differential Evolution

Steered molecular dynamics (SMD) simulation is a powerful method in computer-aided drug design as it can be used to access the relative binding affinity with high precision but with low computational cost. The success of SMD depends on the choice of the direction along which the ligand is pulled from the receptor-binding site. In most simulations, the unidirectional pathway was used, but in some cases, this choice resulted in the ligand colliding with the complex surface of the exit tunnel. To overcome this difficulty, several variants of SMD with multidirectional pulling have been proposed, but they are not completely devoid of disadvantages. Here, we have proposed to determine the direction of pulling with a simple scoring function that minimizes the receptor–ligand interaction, and an optimization algorithm called differential evolution is used for energy minimization. The effectiveness of our protocol was demonstrated by finding expulsion pathways of Huperzine A and camphor from the binding site of Torpedo California acetylcholinesterase and P450cam proteins, respectively, and comparing them with the previous results obtained using memetic sampling and random acceleration molecular dynamics. In addition, by applying this protocol to a set of ligands bound with LSD1 (lysine specific demethylase 1), we obtained a much higher correlation between the work of pulling force and experimental data on the inhibition constant IC50 compared to that obtained using the unidirectional approach based on minimal steric hindrance.

The binding of heparin to spike glycoprotein inhibits SARS-CoV-2 infection by three mechanisms

Heparin, a naturally occurring glycosaminoglycan, has been found to have antiviral activity against SARS-CoV-2, the causative virus of COVID-19. To elucidate the mechanistic basis for the antiviral activity of heparin, we investigated the binding of heparin to the SARS-CoV-2 spike glycoprotein by means of sliding window docking, molecular dynamics simulations, and biochemical assays. Our simulations show that heparin binds at long, positively-charged patches on the spike glycoprotein, thereby masking basic residues of both the receptor binding domain (RBD) and the multifunctional S1/S2 site. Biochemical experiments corroborated the simulation results, showing that heparin inhibits the furin-mediated cleavage of spike by binding to the S1/S2 site. Our simulations also showed that heparin can act on the hinge region responsible for motion of the RBD between the inactive closed and active open conformations of the spike glycoprotein. In simulations of the closed spike homotrimer, heparin binds the RBD and the N-terminal domain of two adjacent spike subunits and hinders opening. In simulations of open spike conformations, heparin induces stabilization of the hinge region and a change in RBD motion. Taken together, our results indicate that heparin can inhibit SARS-CoV-2 infection by three mechanisms: by allosterically hindering binding to the host cell receptor, by directly competing with binding to host heparan sulfate proteoglycan co-receptors, and by preventing spike cleavage by furin. Furthermore, these simulations provide insights into how host heparan sulfate proteoglycans can facilitate viral infection. Our results will aid the rational optimization of heparin derivatives for SARS-CoV-2 antiviral therapy.

The prediction of protein-ligand unbinding for modern drug discovery

INTRODUCTION

Drug-target thermodynamic and kinetic information have perennially important roles in drug design. The prediction of protein-ligand unbinding, which can provide important kinetic information, in experiments continues to face great challenges. Uncovering protein-ligand unbinding through molecular dynamics simulations has become efficient and inexpensive with the progress and enhancement of computing power and sampling methods.

AREAS COVERED

In this review, various sampling methods for protein-ligand unbinding and their basic principles are firstly briefly introduced. Then, their applications in predicting aspects of protein-ligand unbinding, including unbinding pathways, dissociation rate constants, residence time and binding affinity, are discussed.

EXPERT OPINION

Although various sampling methods have been successfully applied in numerous systems, they still have shortcomings and deficiencies. Most enhanced sampling methods require researchers to possess a wealth of prior knowledge of collective variables or reaction coordinates. In addition, most systems studied at present are relatively simple, and the study of complex systems in real drug research remains greatly challenging. Through the combination of machine learning and enhanced sampling methods, prediction accuracy can be further improved, and some problems encountered in complex systems also may be solved.

Protein-Ligand Dissociation Rate Constant from All-Atom Simulation

Dissociation of a ligand isoniazid from a protein catalase was investigated using all-atom molecular dynamics (MD) simulations. Random acceleration MD (τ-RAMD) was used, in which a random artificial force applied to the ligand facilitates its dissociation. We have suggested a novel approach to extrapolate such obtained dissociation times to the zero-force limit assuming never before attempted universal exponential dependence of the bond strength on the applied force, allowing direct comparison with experimentally measured values. We have found that our calculated dissociation time was equal to 36.1 s with statistically significant values distributed in the interval of 0.2-72.0 s, which quantitatively matches the experimental value of 50 ± 8 s despite the extrapolation over 9 orders of magnitude in time.

G Protein-Coupled Receptor-Ligand Dissociation Rates and Mechanisms from τRAMD Simulations

There is a growing appreciation of the importance of drug-target binding kinetics for lead optimization. For G protein-coupled receptors (GPCRs), which mediate signaling over a wide range of time scales, the drug dissociation rate is often a better predictor of in vivo efficacy than binding affinity, although it is more challenging to compute. Here, we assess the ability of the τ-Random Acceleration Molecular Dynamics (τRAMD) approach to reproduce relative residence times and reveal dissociation mechanisms and the effects of allosteric modulation for two important membrane-embedded drug targets: the β2-adrenergic receptor and the muscarinic acetylcholine receptor M2. The dissociation mechanisms observed in the relatively short RAMD simulations (in which molecular dynamics (MD) simulations are performed using an additional force with an adaptively assigned random orientation applied to the ligand) are in general agreement with much more computationally intensive conventional MD and metadynamics simulations. Remarkably, although decreasing the magnitude of the random force generally reduces the number of egress routes observed, the ranking of ligands by dissociation rate is hardly affected and agrees well with experiment. The simulations also reproduce changes in residence time due to allosteric modulation and reveal associated changes in ligand dissociation pathways.

Contact Map Fingerprints of Protein-Ligand Unbinding Trajectories Reveal Mechanisms Determining Residence Times Computed from Scaled Molecular Dynamics

The binding kinetic properties of potential drugs may significantly influence their subsequent clinical efficacy. Predictions of these properties based on computer simulations provide a useful alternative to their expensive and time-consuming experimental counterparts, even at an early drug discovery stage. Herein, we perform scaled molecular dynamics (ScaledMD) simulations on a set of 27 ligands of HSP90 belonging to more than seven chemical series to estimate their relative residence times. We introduce two new techniques for the analysis and the classification of the simulated unbinding trajectories. The first technique, which helps in estimating the limits of the free energy well around the bound state, and the second one, based on a new contact map fingerprint, allow the description and the comparison of the paths that lead to unbinding. Using these analyses, we find that ScaledMD's relative residence time generally enables the identification of the slowest unbinders. We propose an explanation for the underestimation of the residence times of a subset of compounds, and we investigate how the biasing in ScaledMD can affect the mechanistic insights that can be gained from the simulations.

ProLIF: a library to encode molecular interactions as fingerprints

Interaction fingerprints are vector representations that summarize the three-dimensional nature of interactions in molecular complexes, typically formed between a protein and a ligand. This kind of encoding has found many applications in drug-discovery projects, from structure-based virtual-screening to machine-learning. Here, we present ProLIF, a Python library designed to generate interaction fingerprints for molecular complexes extracted from molecular dynamics trajectories, experimental structures, and docking simulations. It can handle complexes formed of any combination of ligand, protein, DNA, or RNA molecules. The available interaction types can be fully reparametrized or extended by user-defined ones. Several tutorials that cover typical use-case scenarios are available, and the documentation is accompanied with code snippets showcasing the integration with other data-analysis libraries for a more seamless user-experience. The library can be freely installed from our GitHub repository (https://github.com/chemosim-lab/ProLIF).

Computer-aided understanding and engineering of enzymatic selectivity

Enzymes offering chemo-, regio-, and stereoselectivity enable the asymmetric synthesis of high-value chiral molecules. Unfortunately, the drawback that naturally occurring enzymes are often inefficient or have undesired selectivity toward non-native substrates hinders the broadening of biocatalytic applications. To match the demands of specific selectivity in asymmetric synthesis, biochemists have implemented various computer-aided strategies in understanding and engineering enzymatic selectivity, diversifying the available repository of artificial enzymes. Here, given that the entire asymmetric catalytic cycle, involving precise interactions within the active pocket and substrate transport in the enzyme channel, could affect the enzymatic efficiency and selectivity, we presented a comprehensive overview of the computer-aided workflow for enzymatic selectivity. This review includes a mechanistic understanding of enzymatic selectivity based on quantum mechanical calculations, rational design of enzymatic selectivity guided by enzyme-substrate interactions, and enzymatic selectivity regulation via enzyme channel engineering. Finally, we discussed the computational paradigm for designing enzyme selectivity in silico to facilitate the advancement of asymmetric biosynthesis.

Computational Enzyme Engineering Pipelines for Optimized Production of Renewable Chemicals

To enable a sustainable supply of chemicals, novel biotechnological solutions are required that replace the reliance on fossil resources. One potential solution is to utilize tailored biosynthetic modules for the metabolic conversion of CO₂ or organic waste to chemicals and fuel by microorganisms. Currently, it is challenging to commercialize biotechnological processes for renewable chemical biomanufacturing because of a lack of highly active and specific biocatalysts. As experimental methods to engineer biocatalysts are time- and cost-intensive, it is important to establish efficient and reliable computational tools that can speed up the identification or optimization of selective, highly active, and stable enzyme variants for utilization in the biotechnological industry. Here, we review and suggest combinations of effective state-of-the-art software and online tools available for computational enzyme engineering pipelines to optimize metabolic pathways for the biosynthesis of renewable chemicals. Using examples relevant for biotechnology, we explain the underlying principles of enzyme engineering and design and illuminate future directions for automated optimization of biocatalysts for the assembly of synthetic metabolic pathways.

The structure of Prp2 bound to RNA and ADP-BeF3- reveals structural features important for RNA unwinding by DEAH-box ATPases

Noncoding intron sequences present in precursor mRNAs need to be removed prior to translation, and they are excised via the spliceosome, a multimegadalton molecular machine composed of numerous protein and RNA components. The DEAH-box ATPase Prp2 plays a crucial role during pre-mRNA splicing as it ensures the catalytic activation of the spliceosome. Despite high structural similarity to other spliceosomal DEAH-box helicases, Prp2 does not seem to function as an RNA helicase, but rather as an RNA-dependent ribonucleoprotein particle-modifying ATPase. Recent crystal structures of the spliceosomal DEAH-box ATPases Prp43 and Prp22, as well as of the related RNA helicase MLE, in complex with RNA have contributed to a better understanding of how RNA binding and processivity might be achieved in this helicase family. In order to shed light onto the divergent manner of function of Prp2, an N-terminally truncated construct of Chaetomium thermophilum Prp2 was crystallized in the presence of ADP-BeF3- and a poly-U12 RNA. The refined structure revealed a virtually identical conformation of the helicase core compared with the ADP-BeF3-- and RNA-bound structure of Prp43, and only a minor shift of the C-terminal domains. However, Prp2 and Prp43 differ in the hook-loop and a loop of the helix-bundle domain, which interacts with the hook-loop and evokes a different RNA conformation immediately after the 3' stack. On replacing these loop residues in Prp43 by the Prp2 sequence, the unwinding activity of Prp43 was abolished. Furthermore, a putative exit tunnel for the γ-phosphate after ATP hydrolysis could be identified in one of the Prp2 structures.

Structure-kinetic relationship reveals the mechanism of selectivity of FAK inhibitors over PYK2

There is increasing evidence of a significant correlation between prolonged drug-target residence time and increased drug efficacy. Here, we report a structural rationale for kinetic selectivity between two closely related kinases: focal adhesion kinase (FAK) and proline-rich tyrosine kinase 2 (PYK2). We found that slowly dissociating FAK inhibitors induce helical structure at the DFG motif of FAK but not PYK2. Binding kinetic data, high-resolution structures and mutagenesis data support the role of hydrophobic interactions of inhibitors with the DFG-helical region, providing a structural rationale for slow dissociation rates from FAK and kinetic selectivity over PYK2. Our experimental data correlate well with computed relative residence times from molecular simulations, supporting a feasible strategy for rationally optimizing ligand residence times. We suggest that the interplay between the protein structural mobility and ligand-induced effects is a key regulator of the kinetic selectivity of inhibitors of FAK versus PYK2.