Sampling of Organic Solutes in Aqueous and Heterogeneous Environments Using Oscillating Excess Chemical Potentials in Grand Canonical-like Monte Carlo-Molecular Dynamics Simulations

High-Throughput Ligand Dissociation Kinetics Predictions Using Site Identification by Ligand Competitive Saturation

The dissociation or off rate, koff, of a drug molecule has been shown to be more relevant to efficacy than affinity for selected systems, motivating the development of predictive computational methodologies. These are largely based on enhanced-sampling molecular dynamics (MD) simulations that come at a high computational cost limiting their utility for drug design where a large number of ligands need to be evaluated. To overcome this, presented is a combined physics- and machine learning (ML)-based approach that uses the physics-based site identification by ligand competitive saturation (SILCS) method to enumerate potential ligand dissociation pathways and calculate ligand dissociation free-energy profiles along those pathways. The calculated free-energy profiles along with molecular properties are used as features to train ML models, including tree and neural network approaches, to predict koff values. The protocol is developed and validated using 329 ligands for 13 proteins showing robustness of the ML workflow built upon the SILCS physics-based free-energy profiles. The resulting SILCS-Kinetics workflow offers a highly efficient method to study ligand dissociation kinetics, providing a powerful tool to facilitate drug design including the ability to generate quantitative estimates of atomic and functional groups contributions to ligand dissociation.

Potent and Selective Human Constitutive Androstane Receptor Activator DL5055 Facilitates Cyclophosphamide-Based Chemotherapies

Enhancement of the metabolic conversion of cyclophosphamide (CPA) increases its therapeutic effects. Activation of the human constitutive androstane receptor (hCAR) induces CYP2B6, a key enzyme responsible for CPA bioactivation. Based on our previous hCAR activator DL5016, we designed and synthesized a series of new hCAR activators. Compared to DL5016, three new compounds 6i, 6k (DL5055), and 7e, showed significantly improved activating potency for hCAR. Particularly, DL5055 activates hCAR with an EC50 of 0.35 μM and EMAX of 4.3, and does not activate hPXR and other related nuclear receptors. It induced the expression of CYP2B6 and caused the translocation of hCAR from the cytoplasm to the nucleus in human primary hepatocytes. DL5055 also induces the expression of Cyp2b10 (the mouse analog of human CYP2B6) in hCAR-transgenic mice. In addition, it significantly enhances the efficacy of CPA-based chemotherapy regimen, CHOP, in a coculture system and a mouse xenograft model in vivo.

Investigating the Interaction between Excipients and Monoclonal Antibodies PGT121 and N49P9.6-FR-LS: A Comprehensive Analysis

N49P9.6-FR-LS and PGT121 are promising antibodies with significant therapeutic potential against HIV infection, but they are prone to precipitation at concentrations greater than 12 to 13 mg/mL. This study evaluates the influence of six excipients─arginine, alanine, sucrose, trehalose, methionine, and glutamate─on the biophysical stability of antibodies. We employed a comprehensive approach, combining computational mAb-excipient interaction analysis via the site-identification by ligand competitive saturation (SILCS) method with extensive experimental characterization. Our experimental matrix included viscosity measurements across temperature gradients, particle size distribution, zeta potential, pH value, and solution appearance, alongside a short-term stability product study at 30 °C and 65% relative humidity, with assessments at t0 (initial), t1 (14 days), and t2 (28 days). Results indicated that sucrose, arginine, alanine, and trehalose provided varying degrees of stabilization for both antibodies. Conversely, glutamate destabilized PGT121 but stabilized N49P9.6-FR-LS, while methionine had a negative effect on N49P9.6-FR-LS but a positive one on PGT121. SILCS-Biologics analysis suggested that stabilization by these excipients is linked to their ability to occupy regions involved in self-protein interactions. Debye-Hückel-Henry charge calculations further indicated that neutral excipients like sucrose and trehalose could alter mAb charges by affecting buffer binding, influencing aggregation propensity. These findings offer valuable insights for optimizing antibody formulations, ensuring enhanced product stability and therapeutic efficacy for HIV treatment.

Mapping the distribution and affinities of ligand interaction sites on human serum albumin

Ligands in many instances interact with a protein at multiple sites with a range of affinities. In this study, ligand-protein interaction sites on human serum albumin (HSA) are mapped using the site-identification by ligand competitive saturation (SILCS)-Biologics approach in conjunction with hydrogen-deuterium exchange (HDX)-mass spectrometry (MS) experiments. Ligands studied include known HSA binders, ibuprofen and ketoprofen, and compounds arginine, alanine, sucrose, and trehalose, excipients used in therapeutic formulations of protein-based drugs. In addition, the impact of excipient binding to HSA on its stability is investigated through temperature-ramp stability studies monitoring solution viscosity. For the studied ligands, interactions that correspond to known drug-binding sites (DSs) are identified. These include previously identified ibuprofen and ketoprofen interaction sites as well as additional sites and, in the case of the excipients, the ligands are shown to also bind at previously unidentified interaction sites, termed excipient sites (ESs) with 20 or more sites identified for the studied compounds. HDX-MS titrations were used to determine dissociation constants for a subset of the interaction sites for ibuprofen, ketoprofen, arginine, and sucrose, which exhibited Kd values in the low micromolar to millimolar range in satisfactory agreement with SILCS-Biologics predicted affinities, validating the computational approach to identify both high- and low-affinity interaction sites. The stability studies indicate the excipients offer protection at low excipient/protein ratios up to 66 with destabilization occurring at ratios above 132 with the exception of sucrose at the t0 time point, indicating that the more favorable affinities of sucrose seen in the SILCS-Biologics and HDX-MS analyses contribute to protein stabilization. These results indicate that ligands can bind to large numbers of interaction sites on proteins, with those interactions having implications for the development of formulations for therapeutic proteins.

Detection of Putative Ligand Dissociation Pathways in Proteins Using Site-Identification by Ligand Competitive Saturation

Drug efficacy often correlates better with dissociation kinetics than binding affinity alone. To study binding kinetics computationally, it is necessary to identify all of the possible ligand dissociation pathways. The site identification by ligand competitive saturation (SILCS) method involves the precomputation of a set of maps (FragMaps), which describe the free energy landscapes of typical chemical functionalities in and around a target protein or RNA. In the current work, we present and implement a method to use SILCS to identify ligand dissociation pathways, termed "SILCS-Pathway." The A* pathfinding algorithm is utilized to enumerate ligand dissociation pathways between the ligand binding site and the surrounding bulk solvent environment defined on evenly spaced points around the protein based on a Fibonacci lattice. The cost function for the A* algorithm is calculated using the SILCS exclusion maps and the SILCS grid free energy scores, thereby identifying paths that account for local protein flexibility and potential favorable interactions with the ligand. By traversing all evenly distributed bulk solvent points around the protein, we located all possible dissociation pathways and clustered them to identify general ligand unbinding pathways. The procedure is verified by using proteins studied previously with enhanced sampling molecular dynamics (MD) techniques and is shown to be capable of capturing important ligand dissociation routes in a highly computationally efficient manner. The identified pathways will act as the foundation for determining ligand dissociation kinetics using SILCS free energy profiles, which will be described in a subsequent article.

Modeling Ligand Binding Site Water Networks with Site Identification by Ligand Competitive Saturation: Impact on Ligand Binding Orientations and Relative Binding Affinities

Appropriate treatment of water contributions to protein-ligand interactions is a very challenging problem in the context of adequately determining the number of waters to investigate and undertaking conformational sampling of the ligands, the waters, and the surrounding protein. In the present study, an extension of the Site Identification by Ligand Competitive Saturation-Monte Carlo (SILCS-MC) docking approach is presented that enables the determination of the location of water molecules in the binding pocket and their impact on the predicted ligand binding orientation and affinities. The approach, termed SILCS-WATER, involves MC sampling of the ligand along with explicit water molecules in a binding site followed by selection of a subset of waters within specified energetic and distance cutoffs that contribute to ligand binding and orientation. To allow for convergence of both the water and ligand orientations, SILCS-WATER is based on just the overlap of the ligand and water with the SILCS FragMaps and the interaction energy between the waters and ligand. Results show that the SILCS-WATER methodology can capture important waters and improve ligand binding orientations. For 6 of 10 multiple ligand-protein systems, the method improved relative binding affinity prediction against experimental results, with substantial improvements in five systems, when compared to standard SILCS-MC. Improved reproduction of crystallographic ligand binding orientations is shown to be an indicator of when SILCS-WATER will yield improved binding affinity correlations. The method also identifies waters interacting with ligands that occupy unfavorable locations with respect to the protein whose displacement through the appropriate ligand modifications should improve ligand binding affinity. Results are consistent with the binding affinity being modeled as a ligand-water complex interacting with the protein. The presented approach offers new possibilities in revealing water networks and their contributions to the binding orientation and affinity of a ligand for a protein and is anticipated to be of utility for computer-aided drug design.

Grand canonical Monte Carlo and deep learning assisted enhanced sampling to characterize the distribution of Mg2+ and influence of the Drude polarizable force field on the stability of folded states of the twister ribozyme

Molecular dynamics simulations are crucial for understanding the structural and dynamical behavior of biomolecular systems, including the impact of their environment. However, there is a gap between the time scale of these simulations and that of real-world experiments. To address this problem, various enhanced simulation methods have been developed. In addition, there has been a significant advancement of the force fields used for simulations associated with the explicit treatment of electronic polarizability. In this study, we apply oscillating chemical potential grand canonical Monte Carlo and machine learning methods to determine reaction coordinates combined with metadynamics simulations to explore the role of Mg2+ distribution and electronic polarizability in the context of the classical Drude oscillator polarizable force field on the stability of the twister ribozyme. The introduction of electronic polarizability along with the details of the distribution of Mg2+ significantly stabilizes the simulations with respect to sampling the crystallographic conformation. The introduction of electronic polarizability leads to increased stability over that obtained with the additive CHARMM36 FF reported in a previous study, allowing for a distribution of a wider range of ions to stabilize twister. Specific interactions contributing to stabilization are identified, including both those observed in the crystal structures and additional experimentally unobserved interactions. Interactions of Mg2+ with the bases are indicated to make important contributions to stabilization. Notably, the presence of specific interactions between the Mg2+ ions and bases or the non-bridging phosphate oxygens (NBPOs) leads to enhanced dipole moments of all three moieties. Mg2+-NBPO interactions led to enhanced dipoles of the phosphates but, interestingly, not in all the participating ions. The present results further indicate the importance of electronic polarizability in stabilizing RNA in molecular simulations and the complicated nature of the relationship of Mg2+-RNA interactions with the polarization response of the bases and phosphates.

Isolation of phytoconstituents from an extract of Murraya paniculata with cytotoxicity and antioxidant activities and in silico evaluation of their potential to bind to aldose reductase (AKR1B1)

The study on Murraya paniculata (Orange Jasmine) stem bark extract found it to have antioxidant and cytotoxic proper-ties. The structures of the isolated phytoconstituents were determined using NMR spectroscopy. Compounds were evaluated for their potential to be aldose reductase inhibitors using molecular docking and dynamics (MD) simulations. Phytochemical screening of methanolic crude extract was performed from which different fractions of the extract were screened for antioxidant activity using the DPPH radical scavenging method and cytotoxicity using the brine shrimp lethality bioassay. The aqueous fraction showed strong antioxidant activity as compared to the standard butylated hy-droxytoluene, whereas pet ether, dichloromethane, chloroform and methanolic extract exhibited moderate antioxidant activity. Activities in the DPPH assay ranged from 17 to 63 µg/ml and all fractions showed cytotoxic activity. Five identified phytochemical compounds (1-5) include ergosterol endoperoxide (1), the coumarin derivatives 7-methoxy-8-(3-methylbut-2-enyl)-1-benzopyran-2-one (2) and 5,7-dimethoxy-8-(3-methylbut-2-enyl)-1-benzopyran-2-one (3) and a mixture of β-sitosterol (4), and stigmasterol (5). Among them ergosterol endoperoxide has been isolated from the stem bark of the M. paniculata for the first time. MD simulations of the identified compounds indicated their potential to bind to the aldose reductase (AKR1B1) protein. Predicted binding affinities of the compounds based on the site identification the ligand competitive saturation (SILCS) technology was -15.04, -8.85, -9.83, -11.95, and -11.75 kcal/mol for 1 through 5, respectively. The present results are anticipated to lead to further study of the activities of the five compounds including experimental evaluation of their inter-actions with AKR1B1.

First-in-class mitogen-activated protein kinase (MAPK) p38α: MAPK-activated protein kinase 2 dual signal modulator with anti-inflammatory and endothelial-stabilizing properties

We previously identified a small molecule, UM101, predicted to bind to the substrate-binding groove of p38α mitogen-activated protein kinase (MAPK) near the binding site of its proinflammatory substrate, mitogen-activated protein kinase-activated protein kinase (MK)2. UM101 exhibited anti-inflammatory, endothelial-stabilizing, and lung-protective effects. To overcome its limited aqueous solubility and p38α binding affinity, we designed an analog of UM101, GEn-1124, with improved aqueous solubility, stability, and p38α-binding affinity. Compared with UM101, GEn-1124 has 18-fold greater p38α-binding affinity as measured by surface plasmon resonance, 11-fold greater aqueous solubility, enhanced barrier-stabilizing activity in thrombin-stimulated human pulmonary artery endothelial cells in vitro, and greater lung protection in vivo. GEn-1124 improved survival from 10%-40% in murine acute lung injury induced by combined exposure to intratracheal bacterial endotoxin lipopolysaccharide instillation and febrile-range hyperthermia and from 0% to 50% in a mouse influenza pneumonia model. Gene expression analysis by RNASeq in tumor necrosis factor α-treated human pulmonary artery endothelial cells showed that the gene-modifying effects of GEn-1124 were much more restricted to tumor necrosis factor α-inducible genes than those of the catalytic site p38 inhibitor, SB203580. Gene expression pathway analysis, confocal immunofluorescence analysis of p38α and MK2 subcellular trafficking, and surface plasmon resonance analysis of phosphorylated p38α:MK2 binding affinity supports a novel mechanism of action. GEn-1124 destabilizes the activated p38α:MK2 complex and dissociates nuclear export of MK2 and p38α, thereby promoting intranuclear retention and enhanced intranuclear signaling by phosphorylated p38α and accelerated inactivation of p38-free cytosolic MK2 by unopposed phosphatases. SIGNIFICANCE STATEMENT: We describe a novel analog of our first-in-class small molecule modulator of p38α/MK2 signaling targeted to a pocket near the glutamate-aspartate-containing substrate binding domain of p38α, which destabilizes the p38α:MK2 complex without blocking p38 catalytic activity or ablating downstream signaling. The result is a rebalancing of downstream proinflammatory and anti-inflammatory signaling, yielding anti-inflammatory, endothelial-stabilizing, and lung-protective effects with therapeutic potential in acute respiratory distress syndrome.

Combined Physics- and Machine-Learning-Based Method to Identify Druggable Binding Sites Using SILCS-Hotspots

Identifying druggable binding sites on proteins is an important and challenging problem, particularly for cryptic, allosteric binding sites that may not be obvious from X-ray, cryo-EM, or predicted structures. The Site-Identification by Ligand Competitive Saturation (SILCS) method accounts for the flexibility of the target protein using all-atom molecular simulations that include various small molecule solutes in aqueous solution. During the simulations, the combination of protein flexibility and comprehensive sampling of the water and solute spatial distributions can identify buried binding pockets absent in experimentally determined structures. Previously, we reported a method for leveraging the information in the SILCS sampling to identify binding sites (termed Hotspots) of small mono- or bicyclic compounds, a subset of which coincide with known binding sites of drug-like molecules. Here, we build on that physics-based approach and present a ML model for ranking the Hotspots according to the likelihood they can accommodate drug-like molecules (e.g., molecular weight >200 Da). In the independent validation set, which includes various enzymes and receptors, our model recalls 67% and 89% of experimentally validated ligand binding sites in the top 10 and 20 ranked Hotspots, respectively. Furthermore, we show that the model's output Decision Function is a useful metric to predict binding sites and their potential druggability in new targets. Given the utility the SILCS method for ligand discovery and optimization, the tools presented represent an important advancement in the identification of orthosteric and allosteric binding sites and the discovery of drug-like molecules targeting those sites.

Exploring Druggable Binding Sites on the Class A GPCRs Using the Residue Interaction Network and Site Identification by Ligand Competitive Saturation

G protein-coupled receptors (GPCRs) play a central role in cellular signaling and are linked to many diseases. Accordingly, computational methods to explore potential allosteric sites for this class of proteins to facilitate the identification of potential modulators are needed. Importantly, the availability of rich structural data providing the locations of the orthosteric ligands and allosteric modulators targeting different GPCRs allows for the validation of approaches to identify new allosteric binding sites. Here, we validate the combination of two computational techniques, the residue interaction network (RIN) model and the site identification by ligand competitive saturation (SILCS) method, to predict putative allosteric binding sites of class A GPCRs. RIN analysis identifies hub residues that mediate allosteric signaling within a receptor and have a high capacity to alter receptor dynamics upon ligand binding. The known orthosteric (and allosteric) binding sites of 18 distinct class A GPCRs were successfully predicted by RIN through a dataset of 105 crystal structures (91 ligand-bound, 14 unbound) with up to 77.8% (76.9%) sensitivity, 92.5% (95.3%) specificity, 51.9% (50%) precision, and 86.2% (92.4%) accuracy based on the experimental and theoretical binding site data. Moreover, graph spectral analysis of the residue networks revealed that the proposed sites were located at the interfaces of highly interconnected residue clusters with a high ability to coordinate the functional dynamics. Then, we employed the SILCS-Hotspots method to assess the druggability of the novel sites predicted for 7 distinct class A GPCRs that are critical for a variety of diseases. While the known orthosteric and allosteric binding sites are successfully explored by our approach, numerous putative allosteric sites with the potential to bind drug-like molecules are proposed. The computational approach presented here promises to be a highly effective tool to predict putative allosteric sites of GPCRs to facilitate the design of effective modulators.

Identifying and Assessing Putative Allosteric Sites and Modulators for CXCR4 Predicted through Network Modeling and Site Identification by Ligand Competitive Saturation

The chemokine receptor CXCR4 is a critical target for the treatment of several cancer types and HIV-1 infections. While orthosteric and allosteric modulators have been developed targeting its extracellular or transmembrane regions, the intramembrane region of CXCR4 may also include allosteric binding sites suitable for the development of allosteric drugs. To investigate this, we apply the Gaussian Network Model (GNM) to the monomeric and dimeric forms of CXCR4 to identify residues essential for its local and global motions located in the hinge regions of the protein. Residue interaction network (RIN) analysis suggests hub residues that participate in allosteric communication throughout the receptor. Mutual residues from the network models reside in regions with a high capacity to alter receptor dynamics upon ligand binding. We then investigate the druggability of these potential allosteric regions using the site identification by ligand competitive saturation (SILCS) approach, revealing two putative allosteric sites on the monomer and three on the homodimer. Two screening campaigns with Glide and SILCS-Monte Carlo docking using FDA-approved drugs suggest 20 putative hit compounds including antifungal drugs, anticancer agents, HIV protease inhibitors, and antimalarial drugs. In vitro assays considering mAB 12G5 and CXCL12 demonstrate both positive and negative allosteric activities of these compounds, supporting our computational approach. However, in vivo functional assays based on the recruitment of β-arrestin to CXCR4 do not show significant agonism and antagonism at a single compound concentration. The present computational pipeline brings a new perspective to computer-aided drug design by combining conformational dynamics based on network analysis and cosolvent analysis based on the SILCS technology to identify putative allosteric binding sites using CXCR4 as a showcase.

Efficient Sampling of Cavity Hydration in Proteins with Nonequilibrium Grand Canonical Monte Carlo and Polarizable Force Fields

Prediction of the hydration levels of protein cavities and active sites is important to both mechanistic analysis and ligand design. Due to the unique microscopic environment of these buried water molecules, a polarizable model is expected to be crucial for an accurate treatment of protein internal hydration in simulations. Here we adapt a nonequilibrium candidate Monte Carlo approach for conducting grand canonical Monte Carlo simulations with the Drude polarizable force field. The GPU implementation enables the efficient sampling of internal cavity hydration levels in biomolecular systems. We also develop an enhanced sampling approach referred to as B-walking, which satisfies detailed balance and readily combines with grand canonical integration to efficiently calculate quantitative binding free energies of water to protein cavities. Applications of these developments are illustrated in a solvent box and the polar ligand binding site in trypsin. Our simulation results show that including electronic polarization leads to a modest but clear improvement in the description of water position and occupancy compared to the crystal structure. The B-walking approach enhances the range of water sampling in different chemical potential windows and thus improves the accuracy of water binding free energy calculations.

Structural and thermodynamic framework for PIEZO1 modulation by small molecules

Mechanosensitive PIEZO channels constitute potential pharmacological targets for multiple clinical conditions, spurring the search for potent chemical PIEZO modulators. Among them is Yoda1, a widely used synthetic small molecule PIEZO1 activator discovered through cell-based high-throughput screening. Yoda1 is thought to bind to PIEZO1's mechanosensory arm domain, sandwiched between two transmembrane regions near the channel pore. However, how the binding of Yoda1 to this region promotes channel activation remains elusive. Here, we first demonstrate that cross-linking PIEZO1 repeats A and B with disulfide bridges reduces the effects of Yoda1 in a redox-dependent manner, suggesting that Yoda1 acts by perturbing the contact between these repeats. Using molecular dynamics-based absolute binding free energy simulations, we next show that Yoda1 preferentially occupies a deeper, amphipathic binding site with higher affinity in PIEZO1 open state. Using Yoda1's binding poses in open and closed states, relative binding free energy simulations were conducted in the membrane environment, recapitulating structure-activity relationships of known Yoda1 analogs. Through virtual screening of an 8 million-compound library using computed fragment maps of the Yoda1 binding site, we subsequently identified two chemical scaffolds with agonist activity toward PIEZO1. This study supports a pharmacological model in which Yoda1 activates PIEZO1 by wedging repeats A and B, providing a structural and thermodynamic framework for the rational design of PIEZO1 modulators. Beyond PIEZO channels, the three orthogonal computational approaches employed here represent a promising path toward drug discovery in highly heterogeneous membrane protein systems.

Combining SILCS and Artificial Intelligence for High-Throughput Prediction of the Passive Permeability of Drug Molecules

Membrane permeability of drug molecules plays a significant role in the development of new therapeutic agents. Accordingly, methods to predict the passive permeability of drug candidates during a medicinal chemistry campaign offer the potential to accelerate the drug design process. In this work, we combine the physics-based site identification by ligand competitive saturation (SILCS) method and data-driven artificial intelligence (AI) to create a high-throughput predictive model for the passive permeability of druglike molecules. In this study, we present a comparative analysis of four regression models to predict membrane permeabilities of small druglike molecules; of the tested models, Random Forest was the most predictive yielding an R2 of 0.81 for the independent data set. The input feature vector used to train the developed prediction model includes absolute free energy profiles of ligands through a POPC-cholesterol bilayer based on ligand grid free energy (LGFE) profiles obtained from the SILCS approach. The use of the membrane free energy profiles from SILCS offers information on the physical forces contributing to ligand permeability, while the use of AI yields a more predictive model trained on experimental PAMPA permeability data for a collection of 229 molecules. This combination allows for rapid estimations of ligand permeability at a level of accuracy beyond currently available predictive models while offering insights into the contributions of the functional groups in the ligands to the permeability barrier, thereby offering quantitative information to facilitate rational ligand design.

GPU-specific algorithms for improved solute sampling in grand canonical Monte Carlo simulations

The Grand Canonical Monte Carlo (GCMC) ensemble defined by the excess chemical potential, μex , volume, and temperature, in the context of molecular simulations allows for variations in the number of particles in the system. In practice, GCMC simulations have been widely applied for the sampling of rare gasses and water, but limited in the context of larger molecules. To overcome this limitation, the oscillating μex GCMC method was introduced and shown to be of utility for sampling small solutes, such as formamide, propane, and benzene, as well as for ionic species such as monocations, acetate, and methylammonium. However, the acceptance of GCMC insertions is low, and the method is computationally demanding. In the present study, we improved the sampling efficiency of the GCMC method using known cavity-bias and configurational-bias algorithms in the context of GPU architecture. Specifically, for GCMC simulations of aqueous solution systems, the configurational-bias algorithm was extended by applying system partitioning in conjunction with a random interval extraction algorithm, thereby improving the efficiency in a highly parallel computing environment. The method is parallelized on the GPU using CUDA and OpenCL, allowing for the code to run on both Nvidia and AMD GPUs, respectively. Notably, the method is particularly well suited for GPU computing as the large number of threads allows for simultaneous sampling of a large number of configurations during insertion attempts without additional computational overhead. In addition, the partitioning scheme allows for simultaneous insertion attempts for large systems, offering considerable efficiency. Calculations on the BK Channel, a transporter, including a lipid bilayer with over 760,000 atoms, show a speed up of ~53-fold through the use of system partitioning. The improved algorithm is then combined with an enhanced μex oscillation protocol and shown to be of utility in the context of the site-identification by ligand competitive saturation (SILCS) co-solvent sampling approach as illustrated through application to the protein CDK2.

Nitro-benzylideneoxymorphone, a bifunctional mu and delta opioid receptor ligand with high mu opioid receptor efficacy

Introduction: There is a major societal need for analgesics with less tolerance, dependence, and abuse liability. Preclinical rodent studies suggest that bifunctional ligands with both mu (MOPr) and delta (DOPr) opioid peptide receptor activity may produce analgesia with reduced tolerance and other side effects. This study explores the structure-activity relationships (SAR) of our previously reported MOPr/DOPr lead, benzylideneoxymorphone (BOM) with C7-methylene-substituted analogs. Methods: Analogs were synthesized and tested in vitro for opioid receptor binding and efficacy. One compound, nitro-BOM (NBOM, 12) was evaluated for antinociceptive effects in the warm water tail withdrawal assay in C57BL/6 mice. Acute and chronic antinociception was determined, as was toxicologic effects on chronic administration. Molecular modeling experiments were performed using the Site Identification by Ligand Competitive Saturation (SILCS) method. Results: NBOM was found to be a potent MOPr agonist/DOPr partial agonist that produces high-efficacy antinociception. Antinociceptive tolerance was observed, as was weight loss; this toxicity was only observed with NBOM and not with BOM. Modeling supports the hypothesis that the increased MOPr efficacy of NBOM is due to the substituted benzylidene ring occupying a nonpolar region within the MOPr agonist state. Discussion: Though antinociceptive tolerance and non-specific toxicity was observed on repeated administration, NBOM provides an important new tool for understanding MOPr/DOPr pharmacology.

Influence of Mg2+ Distribution on the Stability of Folded States of the Twister Ribozyme Revealed Using Grand Canonical Monte Carlo and Generative Deep Learning Enhanced Sampling

Metal ions, particularly magnesium ions (Mg2+), play a role in stabilizing the tertiary structures of RNA molecules. Theoretical models and experimental techniques show that metal ions can change RNA dynamics and how it transitions through different stages of folding. However, the specific ways in which metal ions contribute to the formation and stabilization of RNA's tertiary structure are not fully understood at the atomic level. Here, we combined oscillating excess chemical potential Grand Canonical Monte Carlo (GCMC) and metadynamics to bias toward the sampling of unfolded states using reaction coordinates generated by machine learning allowing for examination of Mg2+-RNA interactions that contribute to stabilizing folded states of the pseudoknot found in the Twister ribozyme. GCMC is used to sample diverse ion distributions around the RNA with deep learning applied to iteratively generate system-specific reaction coordinates to maximize conformational sampling during metadynamics simulations. Results from 6 μs simulations performed on 9 individual systems indicate that Mg2+ ions play a crucial role in stabilizing the three-dimensional (3D) structure of the RNA by stabilizing specific interactions of phosphate groups or phosphate groups and bases of neighboring nucleotides. While many phosphates are accessible to interactions with Mg2+, it is observed that multiple, specific interactions are required to sample conformations close to the folded state; coordination of Mg2+ at individual specific sites facilitates sampling of folded conformations though unfolding ultimately occurs. It is only when multiple specific interactions occur, including the presence of specific inner-shell cation interactions linking two nucleotides, that conformations close to the folded state are stable. While many of the identified Mg2+ interactions are observed in the X-ray crystal structure of Twister, the present study suggests two new Mg2+ ion sites in the Twister ribozyme that contribute to stabilization. In addition, specific interactions with Mg2+ are observed that destabilize the local RNA structure, a process that may facilitate the folding of RNA into its correct structure.

Integrated Covalent Drug Design Workflow Using Site Identification by Ligand Competitive Saturation

Covalent drug design is an important component in drug discovery. Traditional drugs interact with their target in a reversible equilibrium, while irreversible covalent drugs increase the drug-target interaction duration by forming a covalent bond with targeted residues and thus may offer a more effective therapeutic approach. To facilitate the design of this class of ligands, computational methods can be used to help identify reactive nucleophilic residues, frequently cysteines, on a target protein for covalent binding, to test various warhead groups for their potential reactivities, and to predict noncovalent contributions to binding that can facilitate drug-target interactions that are important for binding specificity. To further aid covalent drug design, we extended a functional group mapping approach based on explicit solvent all-atom molecular simulations (SILCS: site identification by ligand competitive saturation) that intrinsically considers protein flexibility, functional group, and protein desolvation along with functional group-protein interactions. Through docking of a library of representative warhead fragments using SILCS-Monte Carlo (SILCS-MC), reactive cysteines can be correctly identified for proteins being tested. Furthermore, a machine learning model was trained to quantify the effectiveness of various warhead groups for proteins using metrics from SILCS-MC as well as experimental model compound warhead reactivity data. The ability to rank covalent molecular binders with similar warheads using SILCS ligand grid free energy (LGFE) ranking was also tested for several proteins. Based on these tools, an integrated SILCS-based workflow was developed, named SILCS-Covalent, which can both qualitatively and quantitatively inform covalent drug discovery.

Site Identification by Ligand Competitive Saturation-Biologics Approach for Structure-Based Protein Charge Prediction

Protein-based therapeutics typically require high concentrations of the active protein, which can lead to protein aggregation and high solution viscosity. Such solution behaviors can limit the stability, bioavailability, and manufacturability of protein-based therapeutics and are directly influenced by the charge of a protein. Protein charge is a system property affected by its environment, including the buffer composition, pH, and temperature. Thus, the charge calculated by summing the charges of each residue in a protein, as is commonly done in computational methods, may significantly differ from the effective charge of the protein as these calculations do not account for contributions from bound ions. Here, we present an extension of the structure-based approach termed site identification by ligand competitive saturation-biologics (SILCS-Biologics) to predict the effective charge of proteins. The SILCS-Biologics approach was applied on a range of protein targets in different salt environments for which membrane-confined electrophoresis-determined charges were previously reported. SILCS-Biologics maps the 3D distribution and predicted occupancy of ions, buffer molecules, and excipient molecules bound to the protein surface in a given salt environment. Using this information, the effective charge of the protein is predicted such that the concentrations of the ions and the presence of excipients or buffers are accounted for. Additionally, SILCS-Biologics also produces 3D structures of the binding sites of ions on the proteins, which enable further analyses such as the characterization of protein surface charge distribution and dipole moments in different environments. Notable is the capability of the method to account for competition between salts, excipients, and buffers on the calculated electrostatic properties in different protein formulations. Our study demonstrates the ability of the SILCS-Biologics approach to predict the effective charge of proteins and its applicability in uncovering protein-ion interactions and their contributions to protein solubility and function.

Structure-Based Design of Potent Iminosugar Inhibitors of Endoplasmic Reticulum α-Glucosidase I with Anti-SARS-CoV-2 Activity

Enveloped viruses depend on the host endoplasmic reticulum (ER) quality control (QC) machinery for proper glycoprotein folding. The endoplasmic reticulum quality control (ERQC) enzyme α-glucosidase I (α-GluI) is an attractive target for developing broad-spectrum antivirals. We synthesized 28 inhibitors designed to interact with all four subsites of the α-GluI active site. These inhibitors are derivatives of the iminosugars 1-deoxynojirimycin (1-DNJ) and valiolamine. Crystal structures of ER α-GluI bound to 25 1-DNJ and three valiolamine derivatives revealed the basis for inhibitory potency. We established the structure-activity relationship (SAR) and used the Site Identification by Ligand Competitive Saturation (SILCS) method to develop a model for predicting α-GluI inhibition. We screened the compounds against SARS-CoV-2 in vitro to identify those with greater antiviral activity than the benchmark α-glucosidase inhibitor UV-4. These host-targeting compounds are candidates for investigation in animal models of SARS-CoV-2 and for testing against other viruses that rely on ERQC for correct glycoprotein folding.

In Silico Discovery of Small Molecule Modulators Targeting the Achilles' Heel of SARS-CoV-2 Spike Protein

The spike protein of SARS-CoV-2 has been a promising target for developing vaccines and therapeutics due to its crucial role in the viral entry process. Previously reported cryogenic electron microscopy (cryo-EM) structures have revealed that free fatty acids (FFA) bind with SARS-CoV-2 spike protein, stabilizing its closed conformation and reducing its interaction with the host cell target in vitro. Inspired by these, we utilized a structure-based virtual screening approach against the conserved FFA-binding pocket to identify small molecule modulators of SARS-CoV-2 spike protein, which helped us identify six hits with micromolar binding affinities. Further evaluation of their commercially available and synthesized analogs enabled us to discover a series of compounds with better binding affinities and solubilities. Notably, our identified compounds exhibited similar binding affinities against the spike proteins of the prototypic SARS-CoV-2 and a currently circulating Omicron BA.4 variant. Furthermore, the cryo-EM structure of the compound SPC-14 bound spike revealed that SPC-14 could shift the conformational equilibrium of the spike protein toward the closed conformation, which is human ACE2 (hACE2) inaccessible. Our identified small molecule modulators targeting the conserved FFA-binding pocket could serve as the starting point for the future development of broad-spectrum COVID-19 intervention treatments.

CryoEM structures of anion exchanger 1 capture multiple states of inward- and outward-facing conformations

Anion exchanger 1 (AE1, band 3) is a major membrane protein of red blood cells and plays a key role in acid-base homeostasis, urine acidification, red blood cell shape regulation, and removal of carbon dioxide during respiration. Though structures of the transmembrane domain (TMD) of three SLC4 transporters, including AE1, have been resolved previously in their outward-facing (OF) state, no mammalian SLC4 structure has been reported in the inward-facing (IF) conformation. Here we present the cryoEM structures of full-length bovine AE1 with its TMD captured in both IF and OF conformations. Remarkably, both IF-IF homodimers and IF-OF heterodimers were detected. The IF structures feature downward movement in the core domain with significant unexpected elongation of TM11. Molecular modeling and structure guided mutagenesis confirmed the functional significance of residues involved in TM11 elongation. Our data provide direct evidence for an elevator-like mechanism of ion transport by an SLC4 family member.

In silico identification of a β2-adrenoceptor allosteric site that selectively augments canonical β2AR-Gs signaling and function

Activation of β2-adrenoceptors (β2ARs) causes airway smooth muscle (ASM) relaxation and bronchodilation, and β2AR agonists (β-agonists) are front-line treatments for asthma and other obstructive lung diseases. However, the therapeutic efficacy of β-agonists is limited by agonist-induced β2AR desensitization and noncanonical β2AR signaling involving β-arrestin that is shown to promote asthma pathophysiology. Accordingly, we undertook the identification of an allosteric site on β2AR that could modulate the activity of β-agonists to overcome these limitations. We employed the site identification by ligand competitive saturation (SILCS) computational method to comprehensively map the entire 3D structure of in silico-generated β2AR intermediate conformations and identified a putative allosteric binding site. Subsequent database screening using SILCS identified drug-like molecules with the potential to bind to the site. Experimental assays in HEK293 cells (expressing recombinant wild-type human β2AR) and human ASM cells (expressing endogenous β2AR) identified positive and negative allosteric modulators (PAMs and NAMs) of β2AR as assessed by regulation of β-agonist-stimulation of cyclic AMP generation. PAMs/NAMs had no effect on β-agonist-induced recruitment of β-arrestin to β2AR- or β-agonist-induced loss of cell surface expression in HEK293 cells expressing β2AR. Mutagenesis analysis of β2AR confirmed the SILCS identified site based on mutants of amino acids R131, Y219, and F282. Finally, functional studies revealed augmentation of β-agonist-induced relaxation of contracted human ASM cells and bronchodilation of contracted airways. These findings identify a allosteric binding site on the β2AR, whose activation selectively augments β-agonist-induced Gs signaling, and increases relaxation of ASM cells, the principal therapeutic effect of β-agonists.

Spatial requirements for ITAM signaling in an intracellular natural killer cell model membrane

FcγRIIIa-FcεRIγ complexes, upon stimulation by antibodies, cluster to initiate intracellular signaling and activate natural killer (NK) cells. Intracellular signaling involves Lck phosphorylation of ITAMs of each monomer of a FcεRIγ homodimer in a FcγRIIIa-FcεRIγ complex and subsequent binding of two phosphotyrosines (pY) in tandem by a Syk family kinase. However, how FcR clustering triggers ITAM signaling is not resolved. Molecular modeling and dynamics (MD) simulations are applied to generate ensembles of structures of the FcγRIIIa and FcεRIγ homodimeric cytoplasmic tails of FcγRIIIa-FcεRIγ complexes based on the transmembrane helices and cytoplasmic tails spaced 120, 80, and 50 Å apart to model different extents of clustering. Site-identification by ligand competitive saturation method with Monte Carlo sampling (SILCS-MC) is used to model how Lck could phosphorylate a diversity of ITAM conformations. At 80 Å separation between FcγRIIIa-FcεRIγ complexes, Lck can perform multiple phosphorylations on individual and multiple ITAMs across complexes, including potential sequential phosphorylation events. Syk may then potentially bind the two pYs within a single ITAM in tandem in isolated FcγRIIIa-FcεRIγ complexes, as observed in CD3ε and ζ chains of T cell receptors by the Syk family kinase ZAP-70. In addition, at 50 Å separation between complexes, unique to natural killer cells over T cells, Syk could potentially bind in tandem to pYs in different ITAMs across FcγRIIIa-FcεRIγ complexes. Thus, we predict that an ensemble of spatial orientations of the ITAMS of FcγRIIIa-FcεRIγ complexes that occur upon clustering lead to ITAM phosphorylation by Lck and subsequent Syk activity thereby facilitating downstream signaling.

SILCS-RNA: Toward a Structure-Based Drug Design Approach for Targeting RNAs with Small Molecules

RNA molecules can act as potential drug targets in different diseases, as their dysregulated expression or misfolding can alter various cellular processes. Noncoding RNAs account for ∼70% of the human genome, and these molecules can have complex tertiary structures that present a great opportunity for targeting by small molecules. In the present study, the site identification by ligand competitive saturation (SILCS) computational approach is extended to target RNA, termed SILCS-RNA. Extensions to the method include an enhanced oscillating excess chemical potential protocol for the grand canonical Monte Carlo calculations and individual simulations of the neutral and charged solutes from which the SILCS functional group affinity maps (FragMaps) are calculated for subsequent binding site identification and docking calculations. The method is developed and evaluated against seven RNA targets and their reported small molecule ligands. SILCS-RNA provides a detailed characterization of the functional group affinity pattern in the small molecule binding sites, recapitulating the types of functional groups present in the ligands. The developed method is also shown to be useful for identification of new potential binding sites and identifying ligand moieties that contribute to binding, granular information that can facilitate ligand design. However, limitations in the method are evident including the ability to map the regions of binding sites occupied by ligand phosphate moieties and to fully account for the wide range of conformational heterogeneity in RNA associated with binding of different small molecules, emphasizing inherent challenges associated with applying computer-aided drug design methods to RNA. While limitations are present, the current study indicates how the SILCS-RNA approach may enhance drug discovery efforts targeting RNAs with small molecules.

hERG Blockade Prediction by Combining Site Identification by Ligand Competitive Saturation and Physicochemical Properties

Human ether-a-go-go-related gene (hERG) potassium channel is well-known contributor to drug-induced cardiotoxicity and therefore an extremely important target when performing safety assessments of drug candidates. Ligand-based approaches in connection with quantitative structure active relationships (QSAR) analyses have been developed to predict hERG toxicity. Availability of the recent published cryogenic electron microscopy (cryo-EM) structure for the hERG channel opened the prospect for using structure-based simulation and docking approaches for hERG drug liability predictions. In recent time, the idea of combining structure- and ligand-based approaches for modeling hERG drug liability has gained momentum offering improvements in predictability when compared to ligand-based QSAR practices alone. The present article demonstrates uniting the structure-based SILCS (site-identification by ligand competitive saturation) approach in conjunction with physicochemical properties to develop predictive models for hERG blockade. This combination leads to improved model predictability based on Pearson's R and percent correct (represents rank-ordering of ligands) metric for different validation sets of hERG blockers involving diverse chemical scaffold and wide range of pIC50 values. The inclusion of the SILCS structure-based approach allows determination of the hERG region to which compounds bind and the contribution of different chemical moieties in the compounds to blockade, thereby facilitating the rational ligand design to minimize hERG liability.

Accelerated Ligand-Mapping Molecular Dynamics Simulations for the Detection of Recalcitrant Cryptic Pockets and Occluded Binding Sites

The identification and characterization of binding sites is a critical component of structure-based drug design (SBDD). Probe-based/cosolvent molecular dynamics (MD) methods that allow for protein flexibility have been developed to predict ligand binding sites. However, cryptic pockets that appear only upon ligand binding and occluded binding sites with no access to the solvent pose significant challenges to these methods. Here, we report the development of accelerated ligand-mapping MD (aLMMD), which combines accelerated MD with LMMD, for the detection of these challenging binding sites. The method was validated on five proteins with what we term "recalcitrant" cryptic pockets, which are deeply buried pockets that require extensive movement of the protein backbone to expose, and three proteins with occluded binding sites. In all the cases, aLMMD was able to detect and sample the binding sites. Our results suggest that aLMMD could be used as a general approach for the detection of such elusive binding sites in protein targets, thus providing valuable information for SBDD.

Application of Site-Identification by Ligand Competitive Saturation in Computer-Aided Drug Design

Site Identification by Ligand Competitive Saturation (SILCS) is a molecular simulation approach that uses diverse small solutes in aqueous solution to obtain functional group affinity patterns of a protein or other macromolecule. This involves employing a combined Grand Canonical Monte Carlo (GCMC)-molecular dynamics (MD) method to sample the full 3D space of the protein, including deep binding pockets and interior cavities from which functional group free energy maps (FragMaps) are obtained. The information content in the maps, which include contributions from protein flexibilty and both protein and functional group desolvation contributions, can be used in many aspects of the drug discovery process. These include identification of novel ligand binding pockets, including allosteric sites, pharmacophore modeling, prediction of relative protein-ligand binding affinities for database screening and lead optimization efforts, evaluation of protein-protein interactions as well as in the formulation of biologics-based drugs including monoclonal antibodies. The present article summarizes the various tools developed in the context of the SILCS methodology and their utility in computer-aided drug design (CADD) applications, showing how the SILCS toolset can improve the drug-development process on a number of fronts with respect to both accuracy and throughput representing a new avenue of CADD applications.

Cryo-EM structure of the sodium-driven chloride/bicarbonate exchanger NDCBE

SLC4 transporters play significant roles in pH regulation and cellular sodium transport. The previously solved structures of the outward facing (OF) conformation for AE1 (SLC4A1) and NBCe1 (SLC4A4) transporters revealed an identical overall fold despite their different transport modes (chloride/bicarbonate exchange versus sodium-carbonate cotransport). However, the exact mechanism determining the different transport modes in the SLC4 family remains unknown. In this work, we report the cryo-EM 3.4 Å structure of the OF conformation of NDCBE (SLC4A8), which shares transport properties with both AE1 and NBCe1 by mediating the electroneutral exchange of sodium-carbonate with chloride. This structure features a fully resolved extracellular loop 3 and well-defined densities corresponding to sodium and carbonate ions in the tentative substrate binding pocket. Further, we combine computational modeling with functional studies to unravel the molecular determinants involved in NDCBE and SLC4 transport.

Rapid and accurate estimation of protein-ligand relative binding affinities using site-identification by ligand competitive saturation

Predicting relative protein-ligand binding affinities is a central pillar of lead optimization efforts in structure-based drug design. The site identification by ligand competitive saturation (SILCS) methodology is based on functional group affinity patterns in the form of free energy maps that may be used to compute protein-ligand binding poses and affinities. Presented are results obtained from the SILCS methodology for a set of eight target proteins as reported originally in Wang et al. (J. Am. Chem. Soc., 2015, 137, 2695-2703) using free energy perturbation (FEP) methods in conjunction with enhanced sampling and cycle closure corrections. These eight targets have been subsequently studied by many other authors to compare the efficacy of their method while comparing with the outcomes of Wang et al. In this work, we present results for a total of 407 ligands on the eight targets and include specific analysis on the subset of 199 ligands considered previously. Using the SILCS methodology we can achieve an average accuracy of up to 77% and 74% when considering the eight targets with their 199 and 407 ligands, respectively, for rank-ordering ligand affinities as calculated by the percent correct metric. This accuracy increases to 82% and 80%, respectively, when the SILCS atomic free energy contributions are optimized using a Bayesian Markov-chain Monte Carlo approach. We also report other metrics including Pearson's correlation coefficient, Pearlman's predictive index, mean unsigned error, and root mean square error for both sets of ligands. The results obtained for the 199 ligands are compared with the outcomes of Wang et al. and other published works. Overall, the SILCS methodology yields similar or better-quality predictions without a priori need for known ligand orientations in terms of the different metrics when compared to current FEP approaches with significant computational savings while additionally offering quantitative estimates of individual atomic contributions to binding free energies. These results further validate the SILCS methodology as an accurate, computationally efficient tool to support lead optimization and drug discovery.

Identification of multiple substrate binding sites in SLC4 transporters in the outward-facing conformation: Insights into the transport mechanism

Solute carrier family 4 (SLC4) transporters mediate the transmembrane transport of HCO₃^-, CO₃^2-, and Cl^- necessary for pH regulation, transepithelial H⁺/base transport, and ion homeostasis. Substrate transport with varying stoichiometry and specificity is achieved through an exchange mechanism and/or through coupling of the uptake of anionic substrates to typically co-transported Na⁺. Recently solved outward-facing structures of two SLC4 members (human anion exchanger 1 [hAE1] and human electrogenic sodium bicarbonate cotransporter 1 [hNBCe1]) with different transport modes (Cl^-/HCO₃^- exchange versus Na⁺-CO₃^2- symport) revealed highly conserved three-dimensional organization of their transmembrane domains. However, the exact location of the ion binding sites and their protein-ion coordination motifs are still unclear. In the present work, we combined site identification by ligand competitive saturation mapping and extensive molecular dynamics sampling with functional mutagenesis studies which led to the identification of two substrate binding sites (entry and central) in the outward-facing states of hAE1 and hNBCe1. Mutation of residues in the identified binding sites led to impaired transport in both proteins. We also showed that R730 in hAE1 is crucial for anion binding in both entry and central sites, whereas in hNBCe1, a Na⁺ acts as an anchor for CO₃^2- binding to the central site. Additionally, protonation of the central acidic residues (E681 in hAE1 and D754 in hNBCe1) alters the ion dynamics in the permeation cavity and may contribute to the transport mode differences in SLC4 proteins. These results provide a basis for understanding the functional differences between hAE1 and hNBCe1 and may facilitate potential drug development for diseases such as proximal and distal renal tubular acidosis.

Functional Group Distributions, Partition Coefficients, and Resistance Factors in Lipid Bilayers Using Site Identification by Ligand Competitive Saturation

Small molecules such as metabolites and drugs must pass through the membrane of the cell, a barrier primarily comprising phospholipid bilayers and embedded proteins. To better understand the process of passive diffusion, knowledge of the ability of various functional groups to partition across bilayers and the associated energetics would be of utility. In the present study, the site identification by ligand competitive saturation (SILCS) methodology has been applied to sample the distributions of a diverse set of chemical solutes representing the functional groups of small molecules across phospholipid bilayers composed of 0.9:0.1 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine/cholesterol and a mixture of 0.52:0.18:0.3 1,2-dioleoyl-sn-glycero-3-phospho-l-serine/1,2-dioleoyl-sn-glycero-3-phosphocholine/cholesterol used in parallel artificial membrane permeability assay experiments. A combination of oscillating chemical potential grand canonical Monte Carlo and molecular dynamics in the SILCS simulations was applied to achieve solute sampling through the bilayers and surrounding aqueous environment from which the distribution of solutes and the functional groups they represent were obtained. Results show differential distribution of aliphatic versus aromatic groups with the former having increased sampling in the center of the bilayers versus in the region of the glycerol linker for the latter. Variations in the distribution of different polar groups are evident, with large differences between negative acetate and positive methylammonium with accumulation of the polar-neutral and acetate solutes above the bilayer head groups. Conversion of the distributions to absolute free energies allows for a detailed understanding of energetics of functional groups in different regions of the bilayers and for calculation of absolute free-energy profiles of multifunctional drug-like molecules across the bilayers from which partition coefficients and resistance factors suitable for insertion into the homogenous solubility-diffusion equation for calculation of permeability were obtained. Comparisons of the calculated bilayer/solution partition coefficients with 1-octanol/water experimental data for both drug-like molecules and the solutes show overall good agreement, validating the calculated distributions and associated absolute free-energy profiles.

Insights into Glucose-6-phosphate Allosteric Activation of β-Glucosidase A

Second-generation ethanol production involves the use of agricultural and forestry waste as feedstock, being an alternative to the first-generation technology as it relies on low-cost abundant residues and does not affect food agriculture. However, the success of second-generation biorefineries relies on energetically efficient processes and effective enzyme cocktails to convert cellulose into fermentable sugars. β-glucosidases catalyze the last step on the enzymatic hydrolysis of cellulose; however, they are often inhibited by glucose. Previous studies demonstrated that glucose-6-phosphate (G6P) is a positive allosteric modulator of Bacillus polymyxa β-glucosidase A, improving enzymatic efficiency, providing thermoresistance, and imparting glucose tolerance. However, the precise molecular details of G6P-β-glucosidase A interactions have not yet been described so far. We investigated the molecular details of G6P binding into B. polymyxa β-glucosidase A through in silico docking using the site identification by ligand competitive saturation technology followed by site-directed mutagenesis studies, from which an allosteric binding site for G6P was identified. In addition, a mechanistic shift toward the transglycosylation reaction as opposed to hydrolysis was observed in the presence of G6P, suggesting a new role of G6P allosteric modulation of the catalytic activity of β-glucosidase A.

Small molecules inhibitors of the heterogeneous ribonuclear protein A18 (hnRNP A18): a regulator of protein translation and an immune checkpoint

We have identified chemical probes that simultaneously inhibit cancer cell progression and an immune checkpoint. Using the computational Site Identification by Ligand Competitive Saturation (SILCS) technology, structural biology and cell-based assays, we identify small molecules that directly and selectively bind to the RNA Recognition Motif (RRM) of hnRNP A18, a regulator of protein translation in cancer cells. hnRNP A18 recognizes a specific RNA signature motif in the 3′UTR of transcripts associated with cancer cell progression (Trx, VEGF, RPA) and, as shown here, a tumor immune checkpoint (CTLA-4). Post-transcriptional regulation of immune checkpoints is a potential therapeutic strategy that remains to be exploited. The probes target hnRNP A18 RRM in vitro and in cells as evaluated by cellular target engagement. As single agents, the probes specifically disrupt hnRNP A18–RNA interactions, downregulate Trx and CTLA-4 protein levels and inhibit proliferation of several cancer cell lines without affecting the viability of normal epithelial cells. These first-in-class chemical probes will greatly facilitate the elucidation of the underexplored biological function of RNA Binding Proteins (RBPs) in cancer cells, including their effects on proliferation and immune checkpoint activation.

Enhancing water sampling of buried binding sites using nonequilibrium candidate Monte Carlo

Water molecules can be found interacting with the surface and within cavities in proteins. However, water exchange between bulk and buried hydration sites can be slow compared to simulation timescales, thus leading to the inefficient sampling of the locations of water. This can pose problems for free energy calculations for computer-aided drug design. Here, we apply a hybrid method that combines nonequilibrium candidate Monte Carlo (NCMC) simulations and molecular dynamics (MD) to enhance sampling of water in specific areas of a system, such as the binding site of a protein. Our approach uses NCMC to gradually remove interactions between a selected water molecule and its environment, then translates the water to a new region, before turning the interactions back on. This approach of gradual removal of interactions, followed by a move and then reintroduction of interactions, allows the environment to relax in response to the proposed water translation, improving acceptance of moves and thereby accelerating water exchange and sampling. We validate this approach on several test systems including the ligand-bound MUP-1 and HSP90 proteins with buried crystallographic waters removed. We show that our BLUES (NCMC/MD) method enhances water sampling relative to normal MD when applied to these systems. Thus, this approach provides a strategy to improve water sampling in molecular simulations which may be useful in practical applications in drug discovery and biomolecular design.

Discovery of beta-lactamase CMY-10 inhibitors for combination therapy against multi-drug resistant Enterobacteriaceae

β-lactam antibiotics are the most widely used antimicrobial agents since the discovery of benzylpenicillin in the 1920s. Unfortunately, these life-saving antibiotics are vulnerable to inactivation by continuously evolving β-lactamase enzymes that are primary resistance determinants in multi-drug resistant pathogens. The current study exploits the strategy of combination therapeutics and aims at identifying novel β-lactamase inhibitors that can inactivate the β-lactamase enzyme of the pathogen while allowing the β-lactam antibiotic to act against its penicillin-binding protein target. Inhibitor discovery applied the Site-Identification by Ligand Competitive Saturation (SILCS) technology to map the functional group requirements of the β-lactamase CMY-10 and generate pharmacophore models of active site. SILCS-MC, Ligand-grid Free Energy (LGFE) analysis and Machine-learning based random-forest (RF) scoring methods were then used to screen and filter a library of 700,000 compounds. From the computational screens 74 compounds were subjected to experimental validation in which β-lactamase activity assay, in vitro susceptibility testing, and Scanning Electron Microscope (SEM) analysis were conducted to explore their antibacterial potential. Eleven compounds were identified as enhancers while 7 compounds were recognized as inhibitors of CMY-10. Of these, compound 11 showed promising activity in β-lactamase activity assay, in vitro susceptibility testing against ATCC strains (E. coli, E. cloacae, E. agglomerans, E. alvei) and MDR clinical isolates (E. cloacae, E. alvei and E. agglomerans), with synergistic assay indicating its potential as a β-lactam enhancer and β-lactamase inhibitor. Structural similarity search against the active compound 11 yielded 28 more compounds. The majority of these compounds also exhibited β-lactamase inhibition potential and antibacterial activity. The non-β-lactam-based β-lactamase inhibitors identified in the current study have the potential to be used in combination therapy with lactam-based antibiotics against MDR clinical isolates that have been found resistant against last-line antibiotics.

Specificity of Molecular Fragments Binding to S100B versus S100A1 as Identified by NMR and Site Identification by Ligand Competitive Saturation (SILCS)

S100B, a biomarker of malignant melanoma, interacts with the p53 protein and diminishes its tumor suppressor function, which makes this S100 family member a promising therapeutic target for treating malignant melanoma. However, it is a challenge to design inhibitors that are specific for S100B in melanoma versus other S100-family members that are important for normal cellular activities. For example, S100A1 is most similar in sequence and structure to S100B, and this S100 protein is important for normal skeletal and cardiac muscle function. Therefore, a combination of NMR and computer aided drug design (CADD) was used to initiate the design of specific S100B inhibitors. Fragment-based screening by NMR, also termed "SAR by NMR," is a well-established method, and was used to examine spectral perturbations in 2D [1H, 15N]-HSQC spectra of Ca2+-bound S100B and Ca2+-bound S100A1, side-by-side, and under identical conditions for comparison. Of the 1000 compounds screened, two were found to be specific for binding Ca2+-bound S100A1 and four were found to be specific for Ca2+-bound S100B, respectively. The NMR spectral perturbations observed in these six data sets were then used to model how each of these small molecule fragments showed specificity for one S100 versus the other using a CADD approach termed Site Identification by Ligand Competitive Saturation (SILCS). In summary, the combination of NMR and computational approaches provided insight into how S100A1 versus S100B bind small molecules specifically, which will enable improved drug design efforts to inhibit elevated S100B in melanoma. Such a fragment-based approach can be used generally to initiate the design of specific inhibitors for other highly homologous drug targets.

Assessing hERG1 Blockade from Bayesian Machine-Learning-Optimized Site Identification by Ligand Competitive Saturation Simulations

Drug-induced cardiotoxicity is a potentially lethal and yet one of the most common side effects with the drugs in clinical use. Most of the drug-induced cardiotoxicity is associated with an off-target pharmacological blockade of K+ currents carried out by the cardiac Human-Ether-a-go-go-Related (hERG1) potassium channel. There is a compulsory preclinical stage safety assessment for the hERG1 blockade for all classes of drugs, which adds substantially to the cost of drug development. The availability of a high-resolution cryogenic electron microscopy (cryo-EM) structure for the channel in its open/depolarized state solved in 2017 enabled the application of molecular modeling for rapid assessment of drug blockade by molecular docking and simulation techniques. More importantly, if successful, in silico methods may allow a path to lead-compound salvaging by mapping out key block determinants. Here, we report the blind application of the site identification by the ligand competitive saturation (SILCS) protocol to map out druggable/regulatory hotspots in the hERG1 channel available for blockers and activators. The SILCS simulations use small solutes representative of common functional groups to sample the chemical space for the entire protein and its environment using all-atom simulations. The resulting chemical maps, FragMaps, explicitly account for receptor flexibility, protein-fragment interactions, and fragment desolvation penalty allowing for rapid ranking of potential ligands as blockers or nonblockers of hERG1. To illustrate the power of the approach, SILCS was applied to a test set of 55 blockers with diverse chemical scaffolds and pIC50 values measured under uniform conditions. The original SILCS model was based on the all-atom modeling of the hERG1 channel in an explicit lipid bilayer and was further augmented with a Bayesian-optimization/machine-learning (BML) stage employing an independent literature-derived training set of 163 molecules. BML approach was used to determine weighting factors for the FragMaps contributions to the scoring function. pIC50 predictions from the combined SILCS/BML approach to the 55 blockers showed a Pearson correlation (PC) coefficient of >0.535 relative to the experimental data. SILCS/BML model was shown to yield substantially improved performance as compared to commonly used rigid and flexible molecular docking methods for a well-established cohort of hERG1 blockers, where no correlation with experimental data was recorded. SILCS/BML results also suggest that a proper weighting of protonation states of common blockers present at physiological pH is essential for accurate predictions of blocker potency. The precalculated and optimized SILCS FragMaps can now be used for the rapid screening of small molecules for their cardiotoxic potential as well as for exploring alternative binding pockets in the hERG1 channel with applications to the rational design of activators.

Computational Characterization of Antibody-Excipient Interactions for Rational Excipient Selection Using the Site Identification by Ligand Competitive Saturation-Biologics Approach

Protein therapeutics typically require a concentrated protein formulation, which can lead to self-association and/or high viscosity due to protein-protein interaction (PPI). Excipients are often added to improve stability, bioavailability, and manufacturability of the protein therapeutics, but the selection of excipients often relies on trial and error. Therefore, understanding the excipient-protein interaction and its effect on non-specific PPI is important for rational selection of formulation development. In this study, we validate a general workflow based on the site identification by ligand competitive saturation (SILCS) technology, termed SILCS-Biologics, that can be applied to protein therapeutics for rational excipient selection. The National Institute of Standards and Technology monoclonal antibody (NISTmAb) reference along with the CNTO607 mAb is used as model antibody proteins to examine PPIs, and NISTmAb was used to further examine excipient-protein interactions, in silico. Metrics from SILCS include the distribution and predicted affinity of excipients, buffer interactions with the NISTmAb Fab, and the relation of the interactions to predicted PPI. Comparison with a range of experimental data showed multiple SILCS metrics to be predictive. Specifically, the number of favorable sites to which an excipient binds and the number of sites to which an excipient binds that are involved in predicted PPIs correlate with the experimentally determined viscosity. In addition, a combination of the number of binding sites and the predicted binding affinity is indicated to be predictive of relative protein stability. Comparison of arginine, trehalose, and sucrose, all of which give the highest viscosity in combination with analysis of B₂₂ and k_D and the SILCS metrics, indicates that higher viscosities are associated with a low number of predicted binding sites, with lower binding affinity of arginine leading to its anomalously high impact on viscosity. The present study indicates the potential for the SILCS-Biologics approach to be of utility in the rational design of excipients during biologics formulation.

Toward Biotherapeutics Formulation Composition Engineering using Site-Identification by Ligand Competitive Saturation (SILCS)

Formulation of protein-based therapeutics employ advanced formulation and analytical technologies for screening various parameters such as buffer, pH, and excipients. At a molecular level, physico-chemical properties of a protein formulation depend on self-interaction between protein molecules, protein-solvent and protein-excipient interactions. This work describes a novel in silico approach, SILCS-Biologics, for structure-based modeling of protein formulations. SILCS Biologics is based on the Site-Identification by Ligand Competitive Saturation (SILCS) technology and enables modeling of interactions among different components of a formulation at an atomistic level while accounting for protein flexibility. It predicts potential hotspot regions on the protein surface for protein-protein and protein-excipient interactions. Here we apply SILCS-Biologics on a Fab domain of a monoclonal antibody (mAbN) to model Fab-Fab interactions and interactions with three amino acid excipients, namely, arginine HCl, proline and lysine HCl. Experiments on 100 mg/ml formulations of mAbN showed that arginine increased, lysine reduced, and proline did not impact viscosity. We use SILCS-Biologics modeling to explore a structure-based hypothesis for the viscosity modulating effect of these excipients. Current efforts are aimed at further validation of this novel computational framework and expanding the scope to model full mAb and other protein therapeutics.

Contributions and competition of Mg2+ and K+ in folding and stabilization of the Twister ribozyme

Native folded and compact intermediate states of RNA typically involve tertiary structures in the presence of divalent ions such as Mg²⁺ in a background of monovalent ions. In a recent study, we have shown how the presence of Mg²⁺ impacts the transition from partially unfolded to folded states through a "push-pull" mechanism where the ion both favors and disfavors the sampling of specific phosphate-phosphate interactions. To further understand the ion atmosphere of RNA in folded and partially folded states results from atomistic umbrella sampling and oscillating chemical potential grand canonical Monte Carlo/molecular dynamics (GCMC/MD) simulations are used to obtain atomic-level details of the distributions of Mg²⁺ and K⁺ ions around Twister RNA. Results show the presence of 100 mM Mg²⁺ to lead to increased charge neutralization over that predicted by counterion condensation theory. Upon going from partially unfolded to folded states, overall charge neutralization increases at all studied ion concentrations that, while associated with an increase in the number of direct ion-phosphate interactions, is fully accounted for by the monovalent K⁺ ions. Furthermore, K⁺ preferentially interacts with purine N7 atoms of helical regions in partially unfolded states, thereby potentially stabilizing the helical regions. Thus, both secondary helical structures and formation of tertiary structures leads to increased counterion condensation, thereby stabilizing those structural features of Twister. Notably, it is shown that K⁺ can act as a surrogate for Mg²⁺ by participating in specific interactions with nonsequential phosphate pairs that occur in the folded state, explaining the ability of Twister to self-cleave at submillimolar Mg²⁺ concentrations.

Accurate Prediction of GPCR Ligand Binding Affinity with Free Energy Perturbation

The computational prediction of relative binding free energies is a crucial goal for drug discovery, and G protein-coupled receptors (GPCRs) are arguably the most important drug target class. However, they present increased complexity to model compared to soluble globular proteins. Despite breakthroughs, experimental X-ray crystal and cryo-EM structures are challenging to attain, meaning computational models of the receptor and ligand binding mode are sometimes necessary. This leads to uncertainty in understanding ligand-protein binding induced changes such as, water positioning and displacement, side chain positioning, hydrogen bond networks, and the overall structure of the hydration shell around the ligand and protein. In other words, the very elements that define structure activity relationships (SARs) and are crucial for accurate binding free energy calculations are typically more uncertain for GPCRs. In this work we use free energy perturbation (FEP) to predict the relative binding free energies for ligands of two different GPCRs. We pinpoint the key aspects for success such as the important role of key water molecules, amino acid ionization states, and the benefit of equilibration with specific ligands. Initial calculations following typical FEP setup and execution protocols delivered no correlation with experiment, but we show how results are improved in a logical and systematic way. This approach gave, in the best cases, a coefficient of determination (R²) compared with experiment in the range of 0.6-0.9 and mean unsigned errors compared to experiment of 0.6-0.7 kcal/mol. We anticipate that our findings will be applicable to other difficult-to-model protein ligand data sets and be of wide interest for the community to continue improving FE binding energy predictions.

Structure-Based Optimization of 10-DEBC Derivatives as Potent and Selective Pim-1 Kinase Inhibitors

Pim-1 kinase has been widely regarded as an attractive target for anticancer drugs. Here, we reported our continued efforts in structure-based optimization of compound 10-DEBC, a previously identified micromolar Pim-1 inhibitor. Guided by the Site Identification by Ligand Competitive Saturation (SILCS) method, we quickly obtained a series of 10-DEBC derivatives with significantly improved activity and selectivity. In particular, compound 26 exhibited an IC₅₀ value of 0.9 nM, as well as 220- and 8-fold selectivity over Pim-2 and Pim-3 kinases, respectively.

Impact of electronic polarizability on protein-functional group interactions

Interactions of proteins with functional groups are key to their biological functions, making it essential that they be accurately modeled. To investigate the impact of the inclusion of explicit treatment of electronic polarizability in force fields on protein-functional group interactions, the additive CHARMM and Drude polarizable force field are compared in the context of the Site-Identification by Ligand Competitive Saturation (SILCS) simulation methodology from which functional group interaction patterns with five proteins for which experimental binding affinities of multiple ligands are available, were obtained. The explicit treatment of polarizability produces significant differences in the functional group interactions in the ligand binding sites including overall enhanced binding of functional groups to the proteins. This is associated with variations of the dipole moments of solutes representative of functional groups in the binding sites relative to aqueous solution with higher dipole moments systematically occurring in the latter, though exceptions occur with positively charged methylammonium. Such variation indicates the complex, heterogeneous nature of the electronic environments of ligand binding sites and emphasizes the inherent limitation of fixed charged, additive force fields for modeling ligand-protein interactions. These effects yield more defined orientation of the functional groups in the binding pockets and a small, but systematic improvement in the ability of the SILCS method to predict the binding orientation and relative affinities of ligands to their target proteins. Overall, these results indicate that the physical model associated with the explicit treatment of polarizability along with the presence of lone pairs in a force field leads to changes in the nature of the interactions of functional groups with proteins versus that occurring with additive force fields, suggesting the utility of polarizable force fields in obtaining a more realistic understanding of protein-ligand interactions.

Mg2+ Impacts the Twister Ribozyme through Push-Pull Stabilization of Nonsequential Phosphate Pairs

RNA molecules perform a variety of biological functions for which the correct three-dimensional structure is essential, including as ribozymes where they catalyze chemical reactions. Metal ions, especially Mg2+, neutralize these negatively charged nucleic acids and specifically stabilize RNA tertiary structures as well as impact the folding landscape of RNAs as they assume their tertiary structures. Specific binding sites of Mg2+ in folded conformations of RNA have been studied extensively; however, the full range of interactions of the ion with compact intermediates and unfolded states of RNA is challenging to investigate, and the atomic details of the mechanism by which the ion facilitates tertiary structure formation is not fully known. Here, umbrella sampling combined with oscillating chemical potential Grand Canonical Monte Carlo/molecular dynamics simulations are used to capture the energetics and atomic-level details of Mg2+-RNA interactions that occur along an unfolding pathway of the Twister ribozyme. The free energy profiles reveal stabilization of partially unfolded states by Mg2+, as observed in unfolding experiments, with this stabilization being due to increased sampling of simultaneous interactions of Mg2+ with two or more nonsequential phosphate groups. Notably, these results indicate a push-pull mechanism in which the Mg2+-RNA interactions actually lead to destabilization of specific nonsequential phosphate-phosphate interactions (i.e., pushed apart), whereas other interactions are stabilized (i.e., pulled together), a balance that stabilizes unfolded states and facilitates the folding of Twister, including the formation of hydrogen bonds associated with the tertiary structure. This study establishes a better understanding of how Mg2+-ion interactions contribute to RNA structural properties and stability.

Ions Everywhere? Mg2+ in the μ-Opioid GPCR and Atomic Details of Their Impact on Function

Filizola and co-workers have applied a combination of long-time molecular dynamics and oscillating chemical potential grand canonical Monte Carlo/molecular dynamics to investigate the distribution of Mg²⁺ and Na⁺ in the μ-opioid receptor and their impact on its function. Results indicate atomic details of potential mechanisms by which Mg²⁺ leads to increased efficacy of opioid analgesics. The presence of information flow between the extracellular loops and the intracellular region of the G-protein-coupled receptors that interacts with G-proteins in the presence of Mg²⁺ may be a phenomenon occurring in other G-protein-coupled receptors and, therefore, potentially of broad impact.

Identification and characterization of fragment binding sites for allosteric ligand design using the site identification by ligand competitive saturation hotspots approach (SILCS-Hotspots)

BACKGROUND

Fragment-based ligand design is used for the development of novel ligands that target macromolecules, most notably proteins. Central to its success is the identification of fragment binding sites that are spatially adjacent such that fragments occupying those sites may be linked to create drug-like ligands. Current experimental and computational approaches that address this problem typically identify only a limited number of sites as well as use a limited number of fragment types.

METHODS

The site-identification by ligand competitive saturation (SILCS) approach is extended to the identification of fragment bindings sites, with the method termed SILCS-Hotspots. The approach involves precomputation of the SILCS FragMaps following which the identification of Hotspots, performed by identifying of all possible fragment binding sites on the full 3D structure of the protein followed by spatial clustering.

RESULTS

The SILCS-Hotspots approach identifies a large number of sites on the target protein, including many sites not accessible in experimental structures due to low binding affinities and binding sites on the protein interior. The identified sites are shown to recapitulate the location of known drug-like molecules in both allosteric and orthosteric binding sites on seven proteins including the androgen receptor, the CDK2 and Erk5 kinases, PTP1B phosphatase and three GPCRs; the β2-adrenergic, GPR40 fatty-acid binding and M2-muscarinic receptors. Analysis indicates the importance of considering all possible fragment binding sites, and not just those accessible to experimental methods, when identifying novel binding sites and performing ligand design versus just considering the most favorable sites. The approach is shown to identify a larger number of known binding sites of drug-like molecules versus the commonly used FTMap and Fpocket methods.

GENERAL SIGNIFICANCE

The present results indicate the potential utility of the SILCS-Hotspots approach for fragment-based rational design of ligands, including allosteric modulators.

Mechanism of μ-Opioid Receptor-Magnesium Interaction and Positive Allosteric Modulation

In the era of opioid abuse epidemics, there is an increased demand for understanding how opioid receptors can be allosterically modulated to guide the development of more effective and safer opioid therapies. Among the modulators of the μ-opioid (MOP) receptor, which is the pharmacological target for the majority of clinically used opioid drugs, are monovalent and divalent cations. Specifically, the monovalent sodium cation (Na⁺) has been known for decades to affect MOP receptor signaling by reducing agonist binding, whereas the divalent magnesium cation (Mg²⁺) has been shown to have the opposite effect, notwithstanding the presence of sodium chloride. Although ultra-high-resolution opioid receptor crystal structures have revealed a specific Na⁺ binding site and molecular dynamics (MD) simulation studies have supported the idea that this monovalent ion reduces agonist binding by stabilizing the receptor inactive state, the putative binding site of Mg²⁺ on the MOP receptor, as well as the molecular determinants responsible for its positive allosteric modulation of the receptor, are unknown. In this work, we carried out tens of microseconds of all-atom MD simulations to investigate the simultaneous binding of Mg²⁺ and Na⁺ cations to inactive and active crystal structures of the MOP receptor embedded in an explicit lipid-water environment and confirmed adequate sampling of Mg²⁺ ion binding with a grand canonical Monte Carlo MD method. Analyses of these simulations shed light on 1) the preferred binding sites of Mg²⁺ on the MOP receptor, 2) details of the competition between Mg²⁺ and Na⁺ cations for specific sites, 3) estimates of binding affinities, and 4) testable hypotheses of the molecular mechanism underlying the positive allosteric modulation of the MOP receptor by the Mg²⁺ cation.

Optimization and Evaluation of Site-Identification by Ligand Competitive Saturation (SILCS) as a Tool for Target-Based Ligand Optimization

Chemical fragment cosolvent sampling techniques have become a versatile tool in ligand-protein binding prediction. Site-identification by ligand competitive saturation (SILCS) is one such method that maps the distribution of chemical fragments on a protein as free energy fields called FragMaps. Ligands are then simulated via Monte Carlo techniques in the field of the FragMaps (SILCS-MC) to predict their binding conformations and relative affinities for the target protein. Application of SILCS-MC using a number of different scoring schemes and MC sampling protocols against multiple protein targets was undertaken to evaluate and optimize the predictive capability of the method. Seven protein targets and 551 ligands with broad chemical variability were used to evaluate and optimize the model to maximize Pearson's correlation coefficient, Pearlman's predictive index, correct relative binding affinity, and root-mean-square error versus the absolute experimental binding affinities. Across the protein-ligand sets, the relative affinities of the ligands were predicted correctly an average of 69% of the time for the highest overall SILCS protocol. Training the FragMap weighting factors using a Bayesian machine learning (ML) algorithm led to an increase to an average 75% relative correct affinity predictions. Furthermore, once the optimal protocol is identified for a specific protein-ligand system average predictabilities of 76% are achieved. The ML algorithm is successful with small training sets of data (30 or more compounds) due to the use of physically correct FragMap weights as priors. Notably, the 76% correct relative prediction rate is similar to or better than free energy perturbation methods that are significantly computationally more expensive than SILCS. The results further support the utility of SILCS as a powerful and computationally accessible tool to support lead optimization and development in drug discovery.