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.

Martini 3 Force Field Parameters for Protein Lipidation Post-Translational Modifications.

Protein lipidations are vital co/post-translational modifications that tether lipid tails to specific protein amino acids, allowing them to anchor to biological membranes, switch their subcellular localization, and modulate association with other proteins. Such lipidations are thus crucial for multiple biological processes including signal transduction, protein trafficking, and membrane localization and are implicated in various diseases as well. Examples of lipid-anchored proteins include the Ras family of proteins that undergo farnesylation; actin and gelsolin that are myristoylated; phospholipase D that is palmitoylated; glycosylphosphatidylinositol-anchored proteins; and others. Here, we develop parameters for cysteine-targeting farnesylation, geranylgeranylation, and palmitoylation, as well as glycine-targeting myristoylation for the latest version of the Martini 3 coarse-grained force field. The parameters are developed using the CHARMM36m all-atom force field parameters as reference. The behavior of the coarse-grained models is consistent with that of the all-atom force field for all lipidations and reproduces key dynamical and structural features of lipid-anchored peptides, such as the solvent-accessible surface area, bilayer penetration depth, and representative conformations of the anchors. The parameters are also validated in simulations of the lipid-anchored peripheral membrane proteins Rheb and Arf1, after comparison with independent all-atom simulations. The parameters, along with mapping schemes for the popular martinize2 tool, are available for download at 10.5281/zenodo.7849262 and also as supporting information.

Membrane Composition and Raf[CRD]-Membrane Attachment Are Driving Forces for K-Ras4B Dimer Stability.

Ras proteins are membrane-anchored GTPases that regulate key cellular signaling networks. It has been recently shown that different anionic lipid types can affect the properties of Ras in terms of dimerization/clustering on the cell membrane. To understand the effects of anionic lipids on key spatiotemporal properties of dimeric K-Ras4B, we perform all-atom molecular dynamics simulations of the dimer K-Ras4B in the presence and absence of Raf[RBD/CRD] effectors on two model anionic lipid membranes: one containing 78% mol DOPC, 20% mol DOPS, and 2% mol PIP2 and another one with enhanced concentration of anionic lipids containing 50% mol DOPC, 40% mol DOPS, and 10% mol PIP2. Analysis of our results unveils the orientational space of dimeric K-Ras4B and shows that the stability of the dimer is enhanced on the membrane containing a high concentration of anionic lipids in the absence of Raf effectors. This enhanced stability is also observed in the presence of Raf[RBD/CRD] effectors although it is not influenced by the concentration of anionic lipids in the membrane, but rather on the ability of Raf[CRD] to anchor to the membrane. We generate dominant K-Ras4B conformations by Markov state modeling and yield the population of states according to the K-Ras4B orientation on the membrane. For the membrane containing anionic lipids, we observe correlations between the diffusion of K-Ras4B and PIP2 and anchoring of anionic lipids to the Raf[CRD] domain. We conclude that the presence of effectors with the Raf[CRD] domain anchoring on the membrane as well as the membrane composition both influence the conformational stability of the K-Ras4B dimer, enabling the preservation of crucial interface interactions.

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.

Evaluating the Efficiency of the Martini Force Field to Study Protein Dimerization in Aqueous and Membrane Environments.

Protein-protein complex assembly is one of the major drivers of biological response. Understanding the mechanisms of protein oligomerization/dimerization would allow one to elucidate how these complexes participate in biological activities and could ultimately lead to new approaches in designing novel therapeutic agents. However, determining the exact association pathways and structures of such complexes remains a challenge. Here, we use parallel tempering metadynamics simulations in the well-tempered ensemble to evaluate the performance of Martini 2.2P and Martini open-beta 3 (Martini 3) force fields in reproducing the structure and energetics of the dimerization process of membrane proteins and proteins in an aqueous solution in reasonable accuracy and throughput. We find that Martini 2.2P systematically overestimates the free energy of association by estimating large barriers in distinct areas, which likely leads to overaggregation when multiple monomers are present. In comparison, the less viscous Martini 3 results in a systematic underestimation of the free energy of association for proteins in solution, while it performs well in describing the association of membrane proteins. In all cases, the near-native dimer complexes are identified as minima in the free energy surface albeit not always as the lowest minima. In the case of Martini 3, we find that the spurious supramolecular protein aggregation present in Martini 2.2P multimer simulations is alleviated and thus this force field may be more suitable for the study of protein oligomerization. We propose that the use of enhanced sampling simulations with a refined coarse-grained force field and appropriately defined collective variables is a robust approach for studying the protein dimerization process, although one should be cautious of the ranking of energy minima.

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.

Rethinking drug design in the artificial intelligence era.

Artificial intelligence (AI) tools are increasingly being applied in drug discovery. While some protagonists point to vast opportunities potentially offered by such tools, others remain sceptical, waiting for a clear impact to be shown in drug discovery projects. The reality is probably somewhere in-between these extremes, yet it is clear that AI is providing new challenges not only for the scientists involved but also for the biopharma industry and its established processes for discovering and developing new medicines. This article presents the views of a diverse group of international experts on the 'grand challenges' in small-molecule drug discovery with AI and the approaches to address them.

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.

A small-molecule inhibitor of C5 complement protein.

The complement pathway is an important part of the immune system, and uncontrolled activation is implicated in many diseases. The human complement component 5 protein (C5) is a validated drug target within the complement pathway, as an anti-C5 antibody (Soliris) is an approved therapy for paroxysmal nocturnal hemoglobinuria. Here, we report the identification, optimization and mechanism of action for the first small-molecule inhibitor of C5 complement protein.

PathwayMap: Molecular Pathway Association with Self-Normalizing Neural Networks.

Drug discovery suffers from high attrition because compounds initially deemed as promising can later show ineffectiveness or toxicity resulting from a poor understanding of their activity profile. In this work, we describe a deep self-normalizing neural network model for the prediction of molecular pathway association and evaluate its performance, showing an AUC ranging from 0.69 to 0.91 on a set of compounds extracted from ChEMBL and from 0.81 to 0.83 on an external data set provided by Novartis. We finally discuss the applicability of the proposed model in the domain of lead discovery. A usable application is available via PlayMolecule.org .

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.

Building New Bridges between In Vitro and In Vivo in Early Drug Discovery: Where Molecular Modeling Meets Systems Biology.

Cellular drug targets exist within networked function-generating systems whose constituent molecular species undergo dynamic interdependent non-equilibrium state transitions in response to specific perturbations (i.e.. inputs). Cellular phenotypic behaviors are manifested through the integrated behaviors of such networks. However, in vitro data are frequently measured and/or interpreted with empirical equilibrium or steady state models (e.g. Hill, Michaelis-Menten, Briggs-Haldane) relevant to isolated target populations. We propose that cells act as analog computers, "solving" sets of coupled "molecular differential equations" (i.e. represented by populations of interacting species)via "integration" of the dynamic state probability distributions among those populations. Disconnects between biochemical and functional/phenotypic assays (cellular/in vivo) may arise with targetcontaining systems that operate far from equilibrium, and/or when coupled contributions (including target-cognate partner binding and drug pharmacokinetics) are neglected in the analysis of biochemical results. The transformation of drug discovery from a trial-and-error endeavor to one based on reliable design criteria depends on improved understanding of the dynamic mechanisms powering cellular function/dysfunction at the systems level. Here, we address the general mechanisms of molecular and cellular function and pharmacological modulation thereof. We outline a first principles theory on the mechanisms by which free energy is stored and transduced into biological function, and by which biological function is modulated by drug-target binding. We propose that cellular function depends on dynamic counter-balanced molecular systems necessitated by the exponential behavior of molecular state transitions under non-equilibrium conditions, including positive versus negative mass action kinetics and solute-induced perturbations to the hydrogen bonds of solvating water versus kT.

"Addition" and "Subtraction": Selectivity Design for Type II Maternal Embryonic Leucine Zipper Kinase Inhibitors.

While adding the structural features that are more favored by on-target activity is the more common strategy in selectivity optimization, the opposite strategy of subtracting the structural features that contribute more to off-target activity can also be very effective. Reported here is our successful effort of improving the kinase selectivity of type II maternal embryonic leucine zipper kinase inhibitors by applying these two complementary approaches together, which clearly demonstrates the powerful synergy between them.

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.

Collaborating to improve the use of free-energy and other quantitative methods in drug discovery.

In May and August, 2016, several pharmaceutical companies convened to discuss and compare experiences with Free Energy Perturbation (FEP). This unusual synchronization of interest was prompted by Schrödinger's FEP+ implementation and offered the opportunity to share fresh studies with FEP and enable broader discussions on the topic. This article summarizes key conclusions of the meetings, including a path forward of actions for this group to aid the accelerated evaluation, application and development of free energy and related quantitative, structure-based design methods.

Toward the Validation of Maternal Embryonic Leucine Zipper Kinase: Discovery, Optimization of Highly Potent and Selective Inhibitors, and Preliminary Biology Insight.

MELK kinase has been implicated in playing an important role in tumorigenesis. Our previous studies suggested that MELK is involved in the regulation of cell cycle and its genetic depletion leads to growth inhibition in a subset of high MELK-expressing basal-like breast cancer cell lines. Herein we describe the discovery and optimization of novel MELK inhibitors 8a and 8b that recapitulate the cellular effects observed by short hairpin ribonucleic acid (shRNA)-mediated MELK knockdown in cellular models. We also discovered a novel fluorine-induced hydrophobic collapse that locked the ligand in its bioactive conformation and led to a 20-fold gain in potency. These novel pharmacological inhibitors achieved high exposure in vivo and were well tolerated, which may allow further in vivo evaluation.

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.

Implications of Dynamic Occupancy, Binding Kinetics, and Channel Gating Kinetics for hERG Blocker Safety Assessment and Mitigation.

Blockade of the hERG potassium channel prolongs the ventricular action potential (AP) and QT interval, and triggers early after depolarizations (EADs) and torsade de pointes (TdP) arrhythmia. Opinions differ as to the causal relationship between hERG blockade and TdP, the relative weighting of other contributing factors, definitive metrics of preclinical proarrhythmicity, and the true safety margin in humans. Here, we have used in silico techniques to characterize the effects of channel gating and binding kinetics on hERG occupancy, and of blockade on the human ventricular AP. Gating effects differ for compounds that are sterically compatible with closed channels (becoming trapped in deactivated channels) versus those that are incompatible with the closed/closing state, and expelled during deactivation. Occupancies of trappable blockers build to equilibrium levels, whereas those of non-trappable blockers build and decay during each AP cycle. Occupancies of ~83% (non-trappable) versus ~63% (trappable) of open/inactive channels caused EADs in our AP simulations. Overall, we conclude that hERG occupancy at therapeutic exposure levels may be tolerated for nontrappable, but not trappable blockers capable of building to the proarrhythmic occupancy level. Furthermore, the widely used Redfern safety index may be biased toward trappable blockers, overestimating the exposure-IC50 separation in nontrappable cases.