Reactivities of acrylamide warheads toward cysteine targets: a QM/ML approach to covalent inhibitor design.

Covalent inhibition offers many advantages over non-covalent inhibition, but covalent warhead reactivity must be carefully balanced to maintain potency while avoiding unwanted side effects. While warhead reactivities are commonly measured with assays, a computational model to predict warhead reactivities could be useful for several aspects of the covalent inhibitor design process. Studies have shown correlations between covalent warhead reactivities and quantum mechanic (QM) properties that describe important aspects of the covalent reaction mechanism. However, the models from these studies are often linear regression equations and can have limitations associated with their usage. Applications of machine learning (ML) models to predict covalent warhead reactivities with QM descriptors are not extensively seen in the literature. This study uses QM descriptors, calculated at different levels of theory, to train ML models to predict reactivities of covalent acrylamide warheads. The QM/ML models are compared with linear regression models built upon the same QM descriptors and with ML models trained on structure-based features like Morgan fingerprints and RDKit descriptors. Experiments show that the QM/ML models outperform the linear regression models and the structure-based ML models, and literature test sets demonstrate the power of the QM/ML models to predict reactivities of unseen acrylamide warhead scaffolds. Ultimately, these QM/ML models are effective, computationally feasible tools that can expedite the design of new covalent inhibitors.

Leveraging Advanced In Silico Techniques in Early Drug Discovery: A Study of Potent Small-Molecule YAP-TEAD PPI Disruptors.

Disruption of the YAP-TEAD protein-protein interaction is an attractive therapeutic strategy in oncology to suppress tumor progression and cancer metastasis. YAP binds to TEAD at a large flat binding interface (∼3500 Å2) devoid of a well-defined druggable pocket, so it has been difficult to design low-molecular-weight compounds to abrogate this protein-protein interaction directly. Recently, work by Furet and coworkers (ChemMedChem 2022, DOI: 10.1002/cmdc.202200303) reported the discovery of the first class of small molecules able to efficiently disrupt the transcriptional activity of TEAD by binding to a specific interaction site of the YAP-TEAD binding interface. Using high-throughput in silico docking, they identified a virtual screening hit from a hot spot derived from their previously rationally designed peptidic inhibitor. Structure-based drug design efforts led to the optimization of the hit compound into a potent lead candidate. Given advances in rapid high-throughput screening and rational approaches to peptidic ligand discovery for challenging targets, we analyzed the pharmacophore features involved in transferring from the peptidic to small-molecule inhibitor that could enable small-molecule discovery for such targets. Here, we show retrospectively that pharmacophore analysis augmented by solvation analysis of molecular dynamics trajectories can guide the designs, while binding free energy calculations provide greater insight into the binding conformation and energetics accompanying the association event. The computed binding free energy estimates agree well with experimental findings and offer useful insight into structural determinants that influence ligand binding to the TEAD interaction surface, even for such a shallow binding site. Taken together, our results demonstrates the utility of advanced in silico methods in structure-based design efforts for difficult-to-drug targets such as the YAP-TEAD transcription factor complex.

Lipid21: Complex Lipid Membrane Simulations with AMBER.

We extend the modular AMBER lipid force field to include anionic lipids, polyunsaturated fatty acid (PUFA) lipids, and sphingomyelin, allowing the simulation of realistic cell membrane lipid compositions, including raft-like domains. Head group torsion parameters are revised, resulting in improved agreement with NMR order parameters, and hydrocarbon chain parameters are updated, providing a better match with phase transition temperature. Extensive validation runs (0.9 µs per lipid type) show good agreement with experimental measurements. Furthermore, the simulation of raft-like bilayers demonstrates the perturbing effect of increasing PUFA concentrations on cholesterol molecules. The force field derivation is consistent with the AMBER philosophy, meaning it can be easily mixed with protein, small molecule, nucleic acid, and carbohydrate force fields.

Combined Free-Energy Calculation and Machine Learning Methods for Understanding Ligand Unbinding Kinetics.

The determination of drug residence times, which define the time an inhibitor is in complex with its target, is a fundamental part of the drug discovery process. Synthesis and experimental measurements of kinetic rate constants are, however, expensive and time consuming. In this work, we aimed to obtain drug residence times computationally. Furthermore, we propose a novel algorithm to identify molecular design objectives based on ligand unbinding kinetics. We designed an enhanced sampling technique to accurately predict the free-energy profiles of the ligand unbinding process, focusing on the free-energy barrier for unbinding. Our method first identifies unbinding paths determining a corresponding set of internal coordinates (ICs) that form contacts between the protein and the ligand; it then iteratively updates these interactions during a series of biased molecular dynamics (MD) simulations to reveal the ICs that are important for the whole of the unbinding process. Subsequently, we performed finite-temperature string simulations to obtain the free-energy barrier for unbinding using the set of ICs as a complex reaction coordinate. Importantly, we also aimed to enable the further design of drugs focusing on improved residence times. To this end, we developed a supervised machine learning (ML) approach with inputs from unbiased "downhill" trajectories initiated near the transition state (TS) ensemble of the string unbinding path. We demonstrate that our ML method can identify key ligand-protein interactions driving the system through the TS. Some of the most important drugs for cancer treatment are kinase inhibitors. One of these kinase targets is cyclin-dependent kinase 2 (CDK2), an appealing target for anticancer drug development. Here, we tested our method using two different CDK2 inhibitors for the potential further development of these compounds. We compared the free-energy barriers obtained from our calculations with those observed in available experimental data. We highlighted important interactions at the distal ends of the ligands that can be targeted for improved residence times. Our method provides a new tool to determine unbinding rates and to identify key structural features of the inhibitors that can be used as starting points for novel design strategies in drug discovery.

Relative Binding Free-Energy Calculations at Lipid-Exposed Sites: Deciphering Hot Spots.

Relative binding free-energy (RBFE) calculations are experiencing resurgence in the computer-aided drug design of novel small molecules due to performance gains allowed by cutting-edge molecular mechanic force fields and computer hardware. Application of RBFE to soluble proteins is becoming a routine, while recent studies outline necessary steps to successfully apply RBFE at the orthosteric site of membrane-embedded G-protein-coupled receptors (GPCRs). In this work, we apply RBFE to a congeneric series of antagonists that bind to a lipid-exposed, extra-helical site of the P2Y1 receptor. We find promising performance of RBFE, such that it may be applied in a predictive manner on drug discovery programs targeting lipid-exposed sites. Further, by the application of the microkinetic model, binding at a lipid-exposed site can be split into (1) membrane partitioning of the drug molecule followed by (2) binding at the extra-helical site. We find that RBFE can be applied to calculate the free energy of each step, allowing the uncoupling of observed binding free energy from the influence of membrane affinity. This protocol may be used to identify binding hot spots at extra-helical sites and guide drug discovery programs toward optimizing intrinsic activity at the target.

Inhibitor binding influences the protonation states of histidines in SARS-CoV-2 main protease.

The main protease (Mpro) of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is an attractive target for antiviral therapeutics. Recently, many high-resolution apo and inhibitor-bound structures of Mpro, a cysteine protease, have been determined, facilitating structure-based drug design. Mpro plays a central role in the viral life cycle by catalyzing the cleavage of SARS-CoV-2 polyproteins. In addition to the catalytic dyad His41-Cys145, Mpro contains multiple histidines including His163, His164, and His172. The protonation states of these histidines and the catalytic nucleophile Cys145 have been debated in previous studies of SARS-CoV Mpro, but have yet to be investigated for SARS-CoV-2. In this work we have used molecular dynamics simulations to determine the structural stability of SARS-CoV-2 Mpro as a function of the protonation assignments for these residues. We simulated both the apo and inhibitor-bound enzyme and found that the conformational stability of the binding site, bound inhibitors, and the hydrogen bond networks of Mpro are highly sensitive to these assignments. Additionally, the two inhibitors studied, the peptidomimetic N3 and an α-ketoamide, display distinct His41/His164 protonation-state-dependent stabilities. While the apo and the N3-bound systems favored N δ (HD) and N ϵ (HE) protonation of His41 and His164, respectively, the α-ketoamide was not stably bound in this state. Our results illustrate the importance of using appropriate histidine protonation states to accurately model the structure and dynamics of SARS-CoV-2 Mpro in both the apo and inhibitor-bound states, a necessary prerequisite for drug-design efforts.

Inhibitor binding influences the protonation states of histidines in SARS-CoV-2 main protease.

The main protease (M pro ) of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is an attractive target for antiviral therapeutics. Recently, many high-resolution apo and inhibitor-bound structures of M pro , a cysteine protease, have been determined, facilitating structure-based drug design. M pro plays a central role in the viral life cycle by catalyzing the cleavage of SARS-CoV-2 polyproteins. In addition to the catalytic dyad His41-Cys145, M pro contains multiple histidines including His163, His164, and His172. The protonation states of these histidines and the catalytic nu-cleophile Cys145 have been debated in previous studies of SARS-CoV M pro , but have yet to be investigated for SARS-CoV-2. In this work we have used molecular dynamics simulations to determine the structural stability of SARS-CoV-2 M pro as a function of the protonation assignments for these residues. We simulated both the apo and inhibitor-bound enzyme and found that the conformational stability of the binding site, bound inhibitors, and the hydrogen bond networks of M pro are highly sensitive to these assignments. Additionally, the two inhibitors studied, the peptidomimetic N3 and an α -ketoamide, display distinct His41/His164 protonation-state-dependent stabilities. While the apo and the N3-bound systems favored N δ (HD) and N ϵ (HE) protonation of His41 and His164, respectively, the α -ketoamide was not stably bound in this state. Our results illustrate the importance of using appropriate histidine protonation states to accurately model the structure and dynamics of SARS-CoV-2 M pro in both the apo and inhibitor-bound states, a necessary prerequisite for drug-design efforts.

Revealing Molecular Determinants of hERG Blocker and Activator Binding.

The Kv11.1 potassium channel, encoded by the human ether-a-go-go-related gene (hERG), plays an essential role in the cardiac action potential. hERG blockade by small molecules can induce "torsade de pointes" arrhythmias and sudden death; as such, it is an important off-target to avoid during drug discovery. Recently, a cryo-EM structure of the open channel state of hERG was reported, opening the door to in silico docking analyses and interpretation of hERG structure-activity relationships, with a view to avoiding blocking activity. Despite this, docking directly to this cryo-EM structure has been reported to yield binding modes that are unable to explain known mutagenesis data. In this work, we use molecular dynamics simulations to sample a range of channel conformations and run ensemble docking campaigns at the known hERG binding site below the selectivity filter, composed of the central cavity and the four deep hydrophobic pockets. We identify a hERG conformational state allowing discrimination of blockers vs nonblockers from docking; furthermore, the binding pocket agrees with mutagenesis data, and blocker binding modes fit the hERG blocker pharmacophore. We then use the same protocol to identify a binding pocket in the hERG channel pore for hERG activators, again agreeing with the reported mutagenesis. Our approach may be useful in drug discovery campaigns to prioritize candidate compounds based on hERG liability via virtual docking screens.

Hidden bias in the DUD-E dataset leads to misleading performance of deep learning in structure-based virtual screening.

Recently much effort has been invested in using convolutional neural network (CNN) models trained on 3D structural images of protein-ligand complexes to distinguish binding from non-binding ligands for virtual screening. However, the dearth of reliable protein-ligand x-ray structures and binding affinity data has required the use of constructed datasets for the training and evaluation of CNN molecular recognition models. Here, we outline various sources of bias in one such widely-used dataset, the Directory of Useful Decoys: Enhanced (DUD-E). We have constructed and performed tests to investigate whether CNN models developed using DUD-E are properly learning the underlying physics of molecular recognition, as intended, or are instead learning biases inherent in the dataset itself. We find that superior enrichment efficiency in CNN models can be attributed to the analogue and decoy bias hidden in the DUD-E dataset rather than successful generalization of the pattern of protein-ligand interactions. Comparing additional deep learning models trained on PDBbind datasets, we found that their enrichment performances using DUD-E are not superior to the performance of the docking program AutoDock Vina. Together, these results suggest that biases that could be present in constructed datasets should be thoroughly evaluated before applying them to machine learning based methodology development.

Using Membrane Partitioning Simulations To Predict Permeability of Forty-Nine Drug-Like Molecules.

A simple descriptor calculated from molecular dynamics simulations of the membrane partitioning event is found to correlate well with experimental measurements of passive membrane permeation from the high-throughput MDCK-LE assay using a data set of 49 drug-like molecules. This descriptor approximates the energy cost of translocation across the hydrophobic membrane core (flip-flop), which for many molecules limits permeability. Performance is found to be superior in comparison to calculated properties such as clogP, clogD, or polar surface area. Furthermore, the atomistic simulations provide a structural understanding of the partitioned drug-membrane complex, facilitating medicinal chemistry optimization of membrane permeability.

Calculating Kinetic Rates and Membrane Permeability from Biased Simulations.

We present a simple approach to calculate the kinetic properties of lipid membrane crossing processes from biased molecular dynamics simulations. We demonstrate that by using biased simulations, one can obtain highly accurate kinetic information with significantly reduced computational time with respect to unbiased simulations. We describe how to conveniently calculate the transition rates to enter, cross, and exit the membrane in terms of the mean first passage times. To obtain free energy barriers and relaxation times from biased simulations only, we constructed Markov models using the dynamic histogram analysis method (DHAM). The permeability coefficients that are calculated from the relaxation times are found to correlate highly with experimentally evaluated values. We show that more generally, certain calculated kinetic properties linked to the crossing of the membrane layer (e.g., barrier height and barrier crossing rates) are good indicators of ordering drugs by permeability. Extending the analysis to a 2D Markov model provides a physical description of the membrane crossing mechanism.

Structure-Kinetic Relationships of Passive Membrane Permeation from Multiscale Modeling.

Passive membrane permeation of small molecules is essential to achieve the required absorption, distribution, metabolism, and excretion (ADME) profiles of drug candidates, in particular intestinal absorption and transport across the blood-brain barrier. Computational investigations of this process typically involve either building QSAR models or performing free energy calculations of the permeation event. Although insightful, these methods rarely bridge the gap between computation and experiment in a quantitative manner, and identifying structural insights to apply toward the design of compounds with improved permeability can be difficult. In this work, we combine molecular dynamics simulations capturing the kinetic steps of permeation at the atomistic level with a dynamic mechanistic model describing permeation at the in vitro level, finding a high level of agreement with experimental permeation measurements. Calculation of the kinetic rate constants determining each step in the permeation event allows derivation of structure-kinetic relationships of permeation. We use these relationships to probe the structural determinants of membrane permeation, finding that the desolvation/loss of hydrogen bonding required to leave the membrane partitioned position controls the membrane flip-flop rate, whereas membrane partitioning determines the rate of leaving the membrane.

Contributions of the membrane dipole potential to the function of voltage-gated cation channels and modulation by small molecule potentiators.

The membrane dipole potential (Ψd) constitutes one of three electrical potentials generated by cell membranes. Ψd arises from the unfavorable parallel alignment of phospholipid and water dipoles, and varies in magnitude both longitudinally and laterally across the bilayer according to membrane composition and phospholipid packing density. In this work, we propose that dynamic counter-balancing between Ψd and the transmembrane potential (ΔΨm) governs the conformational state transitions of voltage-gated ion channels. Ψd consists of 1) static outer, and dynamic inner leaflet components (Ψd(extra) and Ψd(intra), respectively); and 2) a transmembrane component (ΔΨd(inner-outer)), ariing from differences in intra- and extracellular leaflet composition. Ψd(intra), which transitions between high and low energy states (Ψd(intra, high) and Ψd(intra, low)) as a function of channel conformation, is transduced by the pore domain. ΔΨd(inner-outer) is transduced by the voltage-sensing (VS) domain in summation with ΔΨm. Potentiation of voltage-gated ion channels is of interest for the treatment of cardiac, neuronal, and other disorders arising from inherited/acquired ion channel dysfunction. Potentiators are widely believed to alter the rates and voltage-dependencies of channel gating transitions by binding to pockets in the membrane-facing and other regions of ion channel targets. Here, we propose that potentiators alter Ψd(intra) and/or Ψd(extra), thereby increasing or decreasing the energy barriers governing channel gating transitions. We used quantum mechanical and molecular dynamics (MD) simulations to predict the overall Ψd-modulating effects of a series of published positive hERG potentiators partitioned into model DOPC bilayers. Our findings suggest a strong correlation between the magnitude of Ψd-lowering and positive hERG potentiation across the series.

Imaging plasma membrane phase behaviour in live cells using a thiophene-based molecular rotor.

Molecular rotors have emerged as versatile probes of microscopic viscosity in lipid bilayers, although it has proved difficult to find probes that stain both phases equally in phase-separated bilayers. Here, we investigate the use of a membrane-targeting viscosity-sensitive fluorophore based on a thiophene moiety with equal affinity for ordered and disordered lipid domains to probe ordering and viscosity within artificial lipid bilayers and live cell plasma membranes.

Uncoupling the Structure-Activity Relationships of ß2 Adrenergic Receptor Ligands from Membrane Binding.

Ligand binding to membrane proteins may be significantly influenced by the interaction of ligands with the membrane. In particular, the microscopic ligand concentration within the membrane surface solvation layer may exceed that in bulk solvent, resulting in overestimation of the intrinsic protein-ligand binding contribution to the apparent/measured affinity. Using published binding data for a set of small molecules with the ß2 adrenergic receptor, we demonstrate that deconvolution of membrane and protein binding contributions allows for improved structure-activity relationship analysis and structure-based drug design. Molecular dynamics simulations of ligand bound membrane protein complexes were used to validate binding poses, allowing analysis of key interactions and binding site solvation to develop structure-activity relationships of ß2 ligand binding. The resulting relationships are consistent with intrinsic binding affinity (corrected for membrane interaction). The successful structure-based design of ligands targeting membrane proteins may require an assessment of membrane affinity to uncouple protein binding from membrane interactions.

Simulation of lipid bilayer self-assembly using all-atom lipid force fields.

In this manuscript we expand significantly on our earlier communication by investigating the bilayer self-assembly of eight different types of phospholipids in unbiased molecular dynamics (MD) simulations using three widely used all-atom lipid force fields. Irrespective of the underlying force field, the lipids are shown to spontaneously form stable lamellar bilayer structures within 1 microsecond, the majority of which display properties in satisfactory agreement with the experimental data. The lipids self-assemble via the same general mechanism, though at formation rates that differ both between lipid types, force fields and even repeats on the same lipid/force field combination. In addition to zwitterionic phosphatidylcholine (PC) and phosphatidylethanolamine (PE) lipids, anionic phosphatidylserine (PS) and phosphatidylglycerol (PG) lipids are represented. To our knowledge this is the first time bilayer self-assembly of phospholipids with negatively charged head groups is demonstrated in all-atom MD simulations.

Imaging phase separation in model lipid membranes through the use of BODIPY based molecular rotors.

In order to fully understand the dynamics of processes within biological lipid membranes, it is necessary to possess an intimate knowledge of the physical state and ordering of lipids within the membrane. Here we report the use of three molecular rotors based on meso-substituted boron-dipyrrin (BODIPY) in combination with fluorescence lifetime spectroscopy to investigate the viscosity and phase behaviour of model lipid bilayers. In phase-separated giant unilamellar vesicles, we visualise both liquid-ordered (Lo) and liquid-disordered (Ld) phases using fluorescence lifetime imaging microscopy (FLIM), determining their associated viscosity values, and investigate the effect of composition on the viscosity of these phases. Additionally, we use molecular dynamics simulations to investigate the orientation of the BODIPY probes within the bilayer, as well as using molecular dynamics simulations and fluorescence correlation spectroscopy (FCS) to compare diffusion coefficients with those predicted from the fluorescence lifetimes of the probes.

All-atom lipid bilayer self-assembly with the AMBER and CHARMM lipid force fields.

This communication reports the first example of spontaneous lipid bilayer formation in unbiased all-atom molecular dynamics (MD) simulations. Using two different lipid force fields we show simulations started from random mixtures of lipids and water in which four different types of phospholipids self-assemble into organized bilayers in under 1 microsecond.

Lipid14: The Amber Lipid Force Field.

The AMBER lipid force field has been updated to create Lipid14, allowing tensionless simulation of a number of lipid types with the AMBER MD package. The modular nature of this force field allows numerous combinations of head and tail groups to create different lipid types, enabling the easy insertion of new lipid species. The Lennard-Jones and torsion parameters of both the head and tail groups have been revised and updated partial charges calculated. The force field has been validated by simulating bilayers of six different lipid types for a total of 0.5 µs each without applying a surface tension; with favorable comparison to experiment for properties such as area per lipid, volume per lipid, bilayer thickness, NMR order parameters, scattering data, and lipid lateral diffusion. As the derivation of this force field is consistent with the AMBER development philosophy, Lipid14 is compatible with the AMBER protein, nucleic acid, carbohydrate, and small molecule force fields.

Further evaluation of quantum chemical methods for the prediction of non-specific binding of positron emission tomography tracers.

The non-specific binding of candidate positron emission tomography (PET) radiotracers causes resulting PET images to have poor contrast and is a key determinant for the success or failure of imaging drugs. Non-specific binding is thought to arise when radiotracers bind to cell membranes and moieties other than their intended target. Our previous preliminary work has proposed the use of the drug-lipid interaction energy descriptor to predict the level of non-specific binding in vivo using a limited set of ten well known PET radiotracers with kinetic modelling data taken from the literature. This work validates and extends the use of the drug-lipid interaction energy descriptor using a new set of twenty-two candidate PET radiotracers with non-specific binding data recently collected at the same imaging centre with consistent methodology. As with the previous set of radiotracers, a significant correlation is found between the quantum chemical drug-lipid interaction energy and in vivo non-specific binding experimental values. In an effort to speed up the calculation process, several semi-empirical quantum chemical methods were assessed for their ability to reproduce the ab initio results. However no single semi-empirical method was found to consistently reproduce the level of correlation achieved with ab initio quantum chemical methods.