AutoDesigner - Core Design, a De Novo Design Algorithm for Chemical Scaffolds: Application to the Design and Synthesis of Novel Selective Wee1 Inhibitors.

The hit identification stage of a drug discovery program generally involves the design of novel chemical scaffolds with desired biological activity against the target(s) of interest. One common approach is scaffold hopping, which is the manual design of novel scaffolds based on known chemical matter. One major limitation of this approach is narrow chemical space exploration, which can lead to difficulties in maintaining or improving biological activity, selectivity, and favorable property space. Another limitation is the lack of preliminary structure-activity relationship (SAR) data around these designs, which could lead to selecting suboptimal scaffolds to advance lead optimization. To address these limitations, we propose AutoDesigner - Core Design (CoreDesign), a de novo scaffold design algorithm. Our approach is a cloud-integrated, de novo design algorithm for systematically exploring and refining chemical scaffolds against biological targets of interest. The algorithm designs, evaluates, and optimizes a vast range, from millions to billions, of molecules in silico, following defined project parameters encompassing structural novelty, physicochemical attributes, potency, and selectivity using active-learning FEP. To validate CoreDesign in a real-world drug discovery setting, we applied it to the design of novel, potent Wee1 inhibitors with improved selectivity over PLK1. Starting from a single known ligand and receptor structure, CoreDesign rapidly explored over 23 billion molecules to identify 1,342 novel chemical series with a mean of 4 compounds per scaffold. To rapidly analyze this large amount of data and prioritize chemical scaffolds for synthesis, we utilize t-Distributed Stochastic Neighbor Embedding (t-SNE) plots of in silico properties. The chemical space projections allowed us to rapidly identify a structurally novel 5-5 fused core meeting all the hit-identification requirements. Several compounds were synthesized and assayed from the scaffold, displaying good potency against Wee1 and excellent PLK1 selectivity. Our results suggest that CoreDesign can significantly speed up the hit-identification process and increase the probability of success of drug discovery campaigns by allowing teams to bring forward high-quality chemical scaffolds derisked by the availability of preliminary SAR.

FEP Protocol Builder: Optimization of Free Energy Perturbation Protocols Using Active Learning.

Significant improvements have been made in the past decade to methods that rapidly and accurately predict binding affinity through free energy perturbation (FEP) calculations. This has been driven by recent advances in small-molecule force fields and sampling algorithms combined with the availability of low-cost parallel computing. Predictive accuracies of ∼1 kcal mol-1 have been regularly achieved, which are sufficient to drive potency optimization in modern drug discovery campaigns. Despite the robustness of these FEP approaches across multiple target classes, there are invariably target systems that do not display expected performance with default FEP settings. Traditionally, these systems required labor-intensive manual protocol development to arrive at parameter settings that produce a predictive FEP model. Due to the (a) relatively large parameter space to be explored, (b) significant compute requirements, and (c) limited understanding of how combinations of parameters can affect FEP performance, manual FEP protocol optimization can take weeks to months to complete, and often does not involve rigorous train-test set splits, resulting in potential overfitting. These manual FEP protocol development timelines do not coincide with tight drug discovery project timelines, essentially preventing the use of FEP calculations for these target systems. Here, we describe an automated workflow termed FEP Protocol Builder (FEP-PB) to rapidly generate accurate FEP protocols for systems that do not perform well with default settings. FEP-PB uses an active-learning workflow to iteratively search the protocol parameter space to develop accurate FEP protocols. To validate this approach, we applied it to pharmaceutically relevant systems where default FEP settings could not produce predictive models. We demonstrate that FEP-PB can rapidly generate accurate FEP protocols for the previously challenging MCL1 system with limited human intervention. We also apply FEP-PB in a real-world drug discovery setting to generate an accurate FEP protocol for the p97 system. FEP-PB is able to generate a more accurate protocol than the expert user, rapidly validating p97 as amenable to free energy calculations. Additionally, through the active-learning workflow, we are able to gain insight into which parameters are most important for a given system. These results suggest that FEP-PB is a robust tool that can aid in rapidly developing accurate FEP protocols and increasing the number of targets that are amenable to the technology.

AutoDesigner, a De Novo Design Algorithm for Rapidly Exploring Large Chemical Space for Lead Optimization: Application to the Design and Synthesis of d-Amino Acid Oxidase Inhibitors.

The lead optimization stage of a drug discovery program generally involves the design, synthesis, and assaying of hundreds to thousands of compounds. The design phase is usually carried out via traditional medicinal chemistry approaches and/or structure-based drug design (SBDD) when suitable structural information is available. Two of the major limitations of this approach are (1) difficulty in rapidly designing potent molecules that adhere to myriad project criteria, or the multiparameter optimization (MPO) problem, and (2) the relatively small number of molecules explored compared to the vast size of chemical space. To address these limitations, we have developed AutoDesigner, a de novo design algorithm. AutoDesigner employs a cloud-native, multistage search algorithm to carry out successive rounds of chemical space exploration and filtering. Millions to billions of virtual molecules are explored and optimized while adhering to a customizable set of project criteria such as physicochemical properties and potency. Additionally, the algorithm only requires a single ligand with measurable affinity and a putative binding model as a starting point, making it amenable to the early stages of an SBDD project where limited data are available. To assess the effectiveness of AutoDesigner, we applied it to the design of novel inhibitors of d-amino acid oxidase (DAO), a target for the treatment of schizophrenia. AutoDesigner was able to generate and efficiently explore over 1 billion molecules to successfully address a variety of project goals. The compounds generated by AutoDesigner that were synthesized and assayed (1) simultaneously met not only physicochemical criteria, clearance, and central nervous system (CNS) penetration (Kp,uu) cutoffs but also potency thresholds and (2) fully utilize structural data to discover and explore novel interactions and a previously unexplored subpocket in the DAO active site. The reported data demonstrate that AutoDesigner can play a key role in accelerating the discovery of novel, potent chemical matter within the constraints of a given drug discovery lead optimization campaign.

From computer-aided drug discovery to computer-driven drug discovery.

Computational chemistry and structure-based design have traditionally been viewed as a subset of tools that could aid acceleration of the drug discovery process, but were not commonly regarded as a driving force in small molecule drug discovery. In the last decade however, there have been dramatic advances in the field, including (1) development of physics-based computational approaches to accurately predict a broad variety of endpoints from potency to solubility, (2) improvements in artificial intelligence and deep learning methods and (3) dramatic increases in computational power with the advent of GPUs and cloud computing, resulting in the ability to explore and accurately profile vast amounts of drug-like chemical space in silico. There have also been simultaneous advancements in structural biology such as cryogenic electron microscopy (cryo-EM) and computational protein-structure prediction, allowing for access to many more high-resolution 3D structures of novel drug-receptor complexes. The convergence of these breakthroughs has positioned structurally-enabled computational methods to be a driving force behind the discovery of novel small molecule therapeutics. This review will give a broad overview of the synergies in recent advances in the fields of computational chemistry, machine learning and structural biology, in particular in the areas of hit identification, hit-to-lead, and lead optimization.

Combining Cloud-Based Free-Energy Calculations, Synthetically Aware Enumerations, and Goal-Directed Generative Machine Learning for Rapid Large-Scale Chemical Exploration and Optimization.

The hit identification process usually involves the profiling of millions to more recently billions of compounds either via traditional experimental high-throughput screens (HTS) or computational virtual high-throughput screens (vHTS). We have previously demonstrated that, by coupling reaction-based enumeration, active learning, and free energy calculations, a similarly large-scale exploration of chemical space can be extended to the hit-to-lead process. In this work, we augment that approach by coupling large scale enumeration and cloud-based free energy perturbation (FEP) profiling with goal-directed generative machine learning, which results in a higher enrichment of potent ideas compared to large scale enumeration alone, while simultaneously staying within the bounds of predefined drug-like property space. We can achieve this by building the molecular distribution for generative machine learning from the PathFinder rules-based enumeration and optimizing for a weighted sum QSAR-based multiparameter optimization function. We examine the utility of this combined approach by designing potent inhibitors of cyclin-dependent kinase 2 (CDK2) and demonstrate a coupled workflow that can (1) provide a 6.4-fold enrichment improvement in identifying <10 nM compounds over random selection and a 1.5-fold enrichment in identifying <10 nM compounds over our previous method, (2) rapidly explore relevant chemical space outside the bounds of commercial reagents, (3) use generative ML approaches to "learn" the SAR from large scale in silico enumerations and generate novel idea molecules for a flexible receptor site that are both potent and within relevant physicochemical space, and (4) produce over 3 000 000 idea molecules and run 1935 FEP simulations, identifying 69 ideas with a predicted IC50 < 10 nM and 358 ideas with a predicted IC50 < 100 nM. The reported data suggest combining both reaction-based and generative machine learning for ideation results in a higher enrichment of potent compounds over previously described approaches and has the potential to rapidly accelerate the discovery of novel chemical matter within a predefined potency and property space.

Reaction-Based Enumeration, Active Learning, and Free Energy Calculations To Rapidly Explore Synthetically Tractable Chemical Space and Optimize Potency of Cyclin-Dependent Kinase 2 Inhibitors.

The hit-to-lead and lead optimization processes usually involve the design, synthesis, and profiling of thousands of analogs prior to clinical candidate nomination. A hit finding campaign may begin with a virtual screen that explores millions of compounds, if not more. However, this scale of computational profiling is not frequently performed in the hit-to-lead or lead optimization phases of drug discovery. This is likely due to the lack of appropriate computational tools to generate synthetically tractable lead-like compounds in silico, and a lack of computational methods to accurately profile compounds prospectively on a large scale. Recent advances in computational power and methods provide the ability to profile much larger libraries of ligands than previously possible. Herein, we report a new computational technique, referred to as "PathFinder", that uses retrosynthetic analysis followed by combinatorial synthesis to generate novel compounds in synthetically accessible chemical space. In this work, the integration of PathFinder-driven compound generation, cloud-based FEP simulations, and active learning are used to rapidly optimize R-groups, and generate new cores for inhibitors of cyclin-dependent kinase 2 (CDK2). Using this approach, we explored >300 000 ideas, performed >5000 FEP simulations, and identified >100 ligands with a predicted IC50 < 100 nM, including four unique cores. To our knowledge, this is the largest set of FEP calculations disclosed in the literature to date. The rapid turnaround time, and scale of chemical exploration, suggests that this is a useful approach to accelerate the discovery of novel chemical matter in drug discovery campaigns.

Acetyl-CoA carboxylase inhibition by ND-630 reduces hepatic steatosis, improves insulin sensitivity, and modulates dyslipidemia in rats.

Simultaneous inhibition of the acetyl-CoA carboxylase (ACC) isozymes ACC1 and ACC2 results in concomitant inhibition of fatty acid synthesis and stimulation of fatty acid oxidation and may favorably affect the morbidity and mortality associated with obesity, diabetes, and fatty liver disease. Using structure-based drug design, we have identified a series of potent allosteric protein-protein interaction inhibitors, exemplified by ND-630, that interact within the ACC phosphopeptide acceptor and dimerization site to prevent dimerization and inhibit the enzymatic activity of both ACC isozymes, reduce fatty acid synthesis and stimulate fatty acid oxidation in cultured cells and in animals, and exhibit favorable drug-like properties. When administered chronically to rats with diet-induced obesity, ND-630 reduces hepatic steatosis, improves insulin sensitivity, reduces weight gain without affecting food intake, and favorably affects dyslipidemia. When administered chronically to Zucker diabetic fatty rats, ND-630 reduces hepatic steatosis, improves glucose-stimulated insulin secretion, and reduces hemoglobin A1c (0.9% reduction). Together, these data suggest that ACC inhibition by representatives of this series may be useful in treating a variety of metabolic disorders, including metabolic syndrome, type 2 diabetes mellitus, and fatty liver disease.

Docking performance of the glide program as evaluated on the Astex and DUD datasets: a complete set of glide SP results and selected results for a new scoring function integrating WaterMap and glide.

Glide SP mode enrichment results for two preparations of the DUD dataset and native ligand docking RMSDs for two preparations of the Astex dataset are presented. Following a best-practices preparation scheme, an average RMSD of 1.140 Å for native ligand docking with Glide SP is computed. Following the same best-practices preparation scheme for the DUD dataset an average area under the ROC curve (AUC) of 0.80 and average early enrichment via the ROC (0.1 %) metric of 0.12 were observed. 74 and 56 % of the 39 best-practices prepared targets showed AUC over 0.7 and 0.8, respectively. Average AUC was greater than 0.7 for all best-practices protein families demonstrating consistent enrichment performance across a broad range of proteins and ligand chemotypes. In both Astex and DUD datasets, docking performance is significantly improved employing a best-practices preparation scheme over using minimally-prepared structures from the PDB. Enrichment results for WScore, a new scoring function and sampling methodology integrating WaterMap and Glide, are presented for four DUD targets, hivrt, hsp90, cdk2, and fxa. WScore performance in early enrichment is consistently strong and all systems examined show AUC > 0.9 and superior early enrichment to DUD best-practices Glide SP results.

Discovery of a Novel Class of d-Amino Acid Oxidase Inhibitors Using the Schrödinger Computational Platform.

d-Serine is a coagonist of the N-methyl d-aspartate (NMDA) receptor, a key excitatory neurotransmitter receptor. In the brain, d-serine is synthesized from its l-isomer by serine racemase and is metabolized by the D-amino acid oxidase (DAO, DAAO). Many studies have linked decreased d-serine concentration and/or increased DAO expression and enzyme activity to NMDA dysfunction and schizophrenia. Thus, it is feasible to employ DAO inhibitors for the treatment of schizophrenia and other indications. Powered by the Schrödinger computational modeling platform, we initiated a research program to identify novel DAO inhibitors with the best-in-class properties. The program execution leveraged an hDAO FEP+ model to prospectively predict compound potency. A new class of DAO inhibitors with desirable properties has been discovered from this endeavor. Our modeling technology on this program has not only enhanced the efficiency of structure-activity relationship development but also helped to identify a previously unexplored subpocket for further optimization.

2-Aryl benzimidazoles: human SCD1-specific stearoyl coenzyme-A desaturase inhibitors.

A series of potent, benzimidazole-based SCD inhibitors which demonstrate selectivity for the hSCD1 enzyme over the hSCD5 isoform are described. The compounds possess suitable cellular activity and pharmacokinetic properties which render them capable of inhibiting liver SCD activity in a mouse pharmacodynamic assay. These 2-aryl benzimidazoles may serve as valuable tools for studying selective hSCD1-inhibition.

Inhibition of Acetyl-CoA Carboxylase by Phosphorylation or the Inhibitor ND-654 Suppresses Lipogenesis and Hepatocellular Carcinoma.

The incidence of hepatocellular carcinoma (HCC) is rapidly increasing due to the prevalence of obesity and non-alcoholic fatty liver disease, but the molecular triggers that initiate disease development are not fully understood. We demonstrate that mice with targeted loss-of-function point mutations within the AMP-activated protein kinase (AMPK) phosphorylation sites on acetyl-CoA carboxylase 1 (ACC1 Ser79Ala) and ACC2 (ACC2 Ser212Ala) have increased liver de novo lipogenesis (DNL) and liver lesions. The same mutation in ACC1 also increases DNL and proliferation in human liver cancer cells. Consistent with these findings, a novel, liver-specific ACC inhibitor (ND-654) that mimics the effects of ACC phosphorylation inhibits hepatic DNL and the development of HCC, improving survival of tumor-bearing rats when used alone and in combination with the multi-kinase inhibitor sorafenib. These studies highlight the importance of DNL and dysregulation of AMPK-mediated ACC phosphorylation in accelerating HCC and the potential of ACC inhibitors for treatment.

Exploring inhibitor binding at the S' subsites of cathepsin L.

We report a series of noncovalent, reversible inhibitors of cathepsin L that have been designed to explore additional binding interactions with the S' subsites. The design was based on our previously reported crystal structure that suggested the possibility of engineering increased interactions with the S' subsites ( Chowdhury et al. J. Med. Chem. 2002, 45, 5321-5329 ). A representative of these new inhibitors has been co-crystallized with mature cathepsin L, and the structure has been solved and refined at 2.2 A. The inhibitors described in this work extend farther into the S' subsites of cathepsins than any inhibitors reported in the literature thus far. These interactions appear to make use of a S3' subsite that can potentially be exploited for enhanced specificity and/or affinity.

Accelerating drug discovery through tight integration of expert molecular design and predictive scoring.

Modeling protein-ligand interactions has been a central goal of computational chemistry for many years. We here review recent progress toward this goal, and highlight the role free energy calculation methods and computational solvent analysis techniques are now having in drug discovery. We further describe recent use of these methodologies to advance two separate drug discovery programs targeting acetyl-CoA carboxylase and tyrosine kinase 2. These examples suggest that tight integration of sophisticated chemistry teams with state-of-the-art computational methods can dramatically improve the efficiency of small molecule drug discovery.

Molecular surface generation using a variable-radius solvent probe.

Protein-ligand binding occurs through interactions at the molecular surface. Hence, a proper description of this surface is essential to our understanding of the process of molecular recognition. Recent studies have noted the inadequacy of using a fixed 1.4 A solvent probe radius to generate the molecular surface. This assumes that water molecules approach all surface atoms at an equal distance irrespective of polarity, which is not the case. To adequately model the protein-water boundary requires that the solvent probe radius change according to the polarity of its contacting atoms, smaller near polar atoms and larger near apolar atoms. To our knowledge, no method currently exists to generate the molecular surface of a protein in this manner. Using a modification of the marching tetrahedra algorithm, we present a method to generate molecular surfaces using a variable radius solvent probe. The resulting surface lacks many of the unrealistic small crevices in nonpolar regions that are found when utilizing an invariant 1.4 A solvent probe, while maintaining the fine detail of the surface at polar regions. On application of the method on a test set of 20 protein structures taken from the Protein Data Bank (PDB), we also find far fewer empty unsolvated cavities that are present when using only a 1.4 A solvent probe, while the majority of solvated and polar cavities is retained. This suggests that the majority of empty cavities previously observed in protein structures might simply be artifacts of the surfacing method. We also find that the variable probe surface can have significant effects on electrostatic calculations by generating a better tuned description of the protein-water boundary. We also examined the binding interfaces of a diverse set of 55 protein-protein complexes. We find that using a variable probe results in an increase in perceived shape complementarity at these sites compared to using a 1.4 A solvent probe. The molecular volume and surface area are geometric values that determine various important properties for macromolecules, and the altered description afforded by a variable solvent probe molecular surface can have significant implications in protein recognition, energetics, folding, and stability calculations.

Inhibition of acetyl-CoA carboxylase suppresses fatty acid synthesis and tumor growth of non-small-cell lung cancer in preclinical models.

Continuous de novo fatty acid synthesis is a common feature of cancer that is required to meet the biosynthetic demands of a growing tumor. This process is controlled by the rate-limiting enzyme acetyl-CoA carboxylase (ACC), an attractive but traditionally intractable drug target. Here we provide genetic and pharmacological evidence that in preclinical models ACC is required to maintain the de novo fatty acid synthesis needed for growth and viability of non-small-cell lung cancer (NSCLC) cells. We describe the ability of ND-646-an allosteric inhibitor of the ACC enzymes ACC1 and ACC2 that prevents ACC subunit dimerization-to suppress fatty acid synthesis in vitro and in vivo. Chronic ND-646 treatment of xenograft and genetically engineered mouse models of NSCLC inhibited tumor growth. When administered as a single agent or in combination with the standard-of-care drug carboplatin, ND-646 markedly suppressed lung tumor growth in the Kras;Trp53-/- (also known as KRAS p53) and Kras;Stk11-/- (also known as KRAS Lkb1) mouse models of NSCLC. These findings demonstrate that ACC mediates a metabolic liability of NSCLC and that ACC inhibition by ND-646 is detrimental to NSCLC growth, supporting further examination of the use of ACC inhibitors in oncology.

Discovery of MK-8831, A Novel Spiro-Proline Macrocycle as a Pan-Genotypic HCV-NS3/4a Protease Inhibitor.

We have been focused on identifying a structurally different next generation inhibitor to MK-5172 (our Ns3/4a protease inhibitor currently under regulatory review), which would achieve superior pangenotypic activity with acceptable safety and pharmacokinetic profile. These efforts have led to the discovery of a novel class of HCV NS3/4a protease inhibitors containing a unique spirocyclic-proline structural motif. The design strategy involved a molecular-modeling based approach, and the optimization efforts on the series to obtain pan-genotypic coverage with good exposures on oral dosing. One of the key elements in this effort was the spirocyclization of the P2 quinoline group, which rigidified and constrained the binding conformation to provide a novel core. A second focus of the team was also to improve the activity against genotype 3a and the key mutant variants of genotype 1b. The rational application of structural chemistry with molecular modeling guided the design and optimization of the structure-activity relationships have resulted in the identification of the clinical candidate MK-8831 with excellent pan-genotypic activity and safety profile.

Novel Quinoline-Based P2-P4 Macrocyclic Derivatives As Pan-Genotypic HCV NS3/4a Protease Inhibitors.

We have previously reported the discovery of our P2-P4 macrocyclic HCV NS3/4a protease inhibitor MK-5172, which in combination with the NS5a inhibitor MK-8742 recently received a breakthrough therapy designation from the US FDA for treatment of chronic HCV infection. Our goal for the next generation NS3/4a inhibitor was to achieve pan-genotypic activity while retaining the pharmacokinetic profile of MK-5172. One of the areas for follow-up investigation involved replacement of the quinoxaline moiety in MK-5172 with a quinoline and studying the effect of substitution at 4-position of the quinoline. The rationale for this effort was based on molecular modeling, which indicated that such modifications would improve interactions with the S2 subsite, in particular with D79. We wish to report herein the discovery of highly potent inhibitors with pan-genotypic activity and an improved profile over MK-5172, especially against gt-3a and A156 mutants.

Crystal structure of StaL, a glycopeptide antibiotic sulfotransferase from Streptomyces toyocaensis.

Over the past decade, antimicrobial resistance has emerged as a major public health crisis. Glycopeptide antibiotics such as vancomycin and teicoplanin are clinically important for the treatment of Gram-positive bacterial infections. StaL is a 3'-phosphoadenosine 5'-phosphosulfate-dependent sulfotransferase capable of sulfating the cross-linked heptapeptide substrate both in vivo and in vitro, yielding the product A47934, a unique teicoplanin-class glycopeptide antibiotic. The sulfonation reaction catalyzed by StaL constitutes the final step in A47934 biosynthesis. Here we report the crystal structure of StaL and its complex with the cofactor product 3'-phosphoadenosine 5'-phosphate. This is only the second prokaryotic sulfotransferase to be structurally characterized. StaL belongs to the large sulfotransferase family and shows higher similarity to cytosolic sulfotransferases (ST) than to the bacterial ST (Stf0). StaL has a novel dimerization motif, different from any other STs that have been structurally characterized. We have also applied molecular modeling to investigate the binding mode of the unique substrate, desulfo-A47934. Based on the structural analysis and modeling results, a series of residues was mutated and kinetically characterized. In addition to the conserved residues (Lys(12), His(67), and Ser(98)), molecular modeling, fluorescence quenching experiments, and mutagenesis studies identified several other residues essential for substrate binding and/or activity, including Trp(34), His(43), Phe(77), Trp(132), and Glu(205).

Solvated interaction energy (SIE) for scoring protein-ligand binding affinities. 1. Exploring the parameter space.

We present a binding free energy function that consists of force field terms supplemented by solvation terms. We used this function to calibrate the solvation model along with the binding interaction terms in a self-consistent manner. The motivation for this approach was that the solute dielectric-constant dependence of calculated hydration gas-to-water transfer free energies is markedly different from that of binding free energies (J. Comput. Chem. 2003, 24, 954). Hence, we sought to calibrate directly the solvation terms in the context of a binding calculation. The five parameters of the model were systematically scanned to best reproduce the absolute binding free energies for a set of 99 protein-ligand complexes. We obtained a mean unsigned error of 1.29 kcal/mol for the predicted absolute binding affinity in a parameter space that was fairly shallow near the optimum. The lowest errors were obtained with solute dielectric values of Din = 20 or higher and scaling of the intermolecular van der Waals interaction energy by factors ranging from 0.03 to 0.15. The high apparent Din and strong van der Waals scaling may reflect the anticorrelation of the change in solvated potential energy and configurational entropy, that is, enthalpy-entropy compensation in ligand binding (Biophys. J. 2004, 87, 3035-3049). Five variations of preparing the protein-ligand data set were explored in order to examine the effect of energy refinement and the presence of bound water on the calculated results. We find that retaining water in the final protein structure used for calculating the binding free energy is not necessary to obtain good results; that is the continuum solvation model is sufficient. Virtual screening enrichment studies on estrogen receptor and thymidine kinase showed a good ability of the binding free energy function to recover true hits in a collection of decoys.

Coupled atomic charge selectivity for optimal ligand-charge distributions at protein binding sites.

Charge optimization as a tool for both analyzing and enhancing binding electrostatics has become an attractive approach over the past few years. An interesting feature of this method for molecular design is that it provides not only the optimal charge magnitudes, but also the selectivity of a particular atomic center for its optimal charge. The current approach to compute the charge selectivity at a given atomic center of a ligand in a particular binding process is based on the binding-energy cost incurred upon the perturbation of the optimal charge distribution by a unit charge at the given atomic center, while keeping the other atomic partial charges at their optimal values. A limitation of this method is that it does not take into account the possible concerted changes in the other atomic charges that may incur a lower energetic cost than perturbing a single charge. Here, we describe a novel approach for characterizing charge selectivity in a concerted manner, taking into account the coupling between the ligand charge centers in the binding process. We apply this novel charge selectivity measure to the celecoxib molecule, a nonsteroidal anti-inflammatory agent binding to cyclooxygenase 2 (COX2), which has been recently shown to also exhibit cross-reactivity toward carbonic anhydrase II (CAII), to which it binds with nanomolar affinity. The uncoupled and coupled charge selectivity profiles over the atomic centers of the celecoxib ligand, binding independently to COX2 and CAII, are analyzed comparatively and rationalized with respect to available experimental data. Very different charge selectivity profiles are obtained for the uncoupled versus coupled selectivity calculations.