Reliably assessing the electronic structure of cytochrome P450 on today's classical computers and tomorrow's quantum computers.

An accurate assessment of how quantum computers can be used for chemical simulation, especially their potential computational advantages, provides important context on how to deploy these future devices. To perform this assessment reliably, quantum resource estimates must be coupled with classical computations attempting to answer relevant chemical questions and to define the classical algorithms simulation frontier. Herein, we explore the quantum computation and classical computation resources required to assess the electronic structure of cytochrome P450 enzymes (CYPs) and thus define a classical-quantum advantage boundary. This is accomplished by analyzing the convergence of density matrix renormalization group plus n-electron valence state perturbation theory (DMRG+NEVPT2) and coupled-cluster singles doubles with noniterative triples [CCSD(T)] calculations for spin gaps in models of the CYP catalytic cycle that indicate multireference character. The quantum resources required to perform phase estimation using qubitized quantum walks are calculated for the same systems. Compilation into the surface code provides runtime estimates to compare directly to DMRG runtimes and to evaluate potential quantum advantage. Both classical and quantum resource estimates suggest that simulation of CYP models at scales large enough to balance dynamic and multiconfigurational electron correlation has the potential to be a quantum advantage problem and emphasizes the important interplay between classical computations and quantum algorithms development for chemical simulation.

BRD4: quantum mechanical protein-ligand binding free energies using the full-protein DFT-based QM-PBSA method.

Fully quantum mechanical approaches to calculating protein-ligand free energies of binding have the potential to reduce empiricism and explicitly account for all physical interactions responsible for protein-ligand binding. In this study, we show a realistic test of the linear-scaling DFT-based QM-PBSA method to estimate quantum mechanical protein-ligand binding free energies for a set of ligands binding to the pharmaceutical drug-target bromodomain containing protein 4 (BRD4). We show that quantum mechanical QM-PBSA is a significant improvement over traditional MM-PBSA in terms of accuracy against experiment and ligand rank ordering and that the quantum and classical binding energies are converged to a similar degree. We test the interaction entropy and normal mode entropy correction terms to QM- and MM-PBSA.

Covalent inhibitor reactivity prediction by the electrophilicity index-in and out of scope.

Drug discovery is an expensive and time-consuming process. To make this process more efficient quantum chemistry methods can be employed. The electrophilicity index is one property that can be calculated by quantum chemistry methods, and if calculated correctly gives insight into the reactivity of covalent inhibitors. Herein we present the usage of the electrophilicity index on three common warheads, i.e., acrylamides, 2-chloroacetamides, and propargylamides. We thoroughly examine the properties of the electrophilicity index, show which pitfalls should be avoided, and what the requirements to successfully apply the electrophilicity index are.

Protein-ligand free energies of binding from full-protein DFT calculations: convergence and choice of exchange-correlation functional.

The accurate prediction of protein-ligand binding free energies with tractable computational methods has the potential to revolutionize drug discovery. Modeling the protein-ligand interaction at a quantum mechanical level, instead of relying on empirical classical-mechanics methods, is an important step toward this goal. In this study, we explore the QM-PBSA method to calculate the free energies of binding of seven ligands to the T4-lysozyme L99A/M102Q mutant using linear-scaling density functional theory on the whole protein-ligand complex. By leveraging modern high-performance computing we perform over 2900 full-protein (2600 atoms) DFT calculations providing new insights into the convergence, precision and reproducibility of the QM-PBSA method. We find that even at moderate sampling over 50 snapshots, the convergence of QM-PBSA is similar to traditional MM-PBSA and that the DFT-based energy evaluations are very reproducible. We show that in the QM-PBSA framework, the physically-motivated GGA exchange-correlation functional PBE outperforms the more modern, dispersion-including non-local and meta-GGA-nonlocal functionals VV10 and B97M-rV. Different empirical dispersion corrections perform similarly well but the three-body dispersion term, as included in Grimme's D3 dispersion, is significant and improves results slightly. Inclusion of an entropy correction term sampled over less than 25 snapshots is detrimental while an entropy correction sampled over the same 50 or 100 snapshots as the enthalpies improves the accuracy of the QM-PBSA method. As full-protein DFT calculations can now be performed on modest computational resources our study demonstrates that they can be a useful addition to the toolbox of free energy calculations.

BIreactive: A Machine-Learning Model to Estimate Covalent Warhead Reactivity.

In the past decade, the pharmaceutical industry has paid closer attention to covalent drugs. Differently from standard noncovalent drugs, these compounds can exhibit peculiar properties, such as higher potency or longer duration of target inhibition with a potentially lower dosage. These properties are mainly driven by the reactive functional group present in the compound, the so-called warhead that forms a covalent bond with a specific nucleophilic amino-acid on the target. In this work, we report the possibility to combine ab initio activation energies with machine-learning to estimate covalent compound intrinsic reactivity. The idea behind this approach is to have a precise estimation of the transition state barriers, and thus of the compound reactivity, but with the speed of a machine-learning algorithm. We call this method "BIreactive". Here, we demonstrate this approach on acrylamides and 2-chloroacetamides, two warhead classes that possess different reaction mechanisms. In combination with our recently implemented truncation algorithm, we also demonstrate the possibility to use BIreactive not only for fragments but also for lead-like molecules. The generic nature of this approach allows also the extension to several other warheads. The combination of these factors makes BIreactive a valuable tool for the covalent drug discovery process in a pharmaceutical context.

A Fast Ab Initio Predictor Tool for Covalent Reactivity Estimation of Acrylamides.

Thanks to their unique mode of action, covalent drugs represent an exceptional opportunity for drug design. After binding to a biologically relevant target system, covalent compounds form a reversible or irreversible covalent bond with a nucleophilic amino acid. Due to the inherently large binding energy of a covalent bond, covalent binders exhibit higher potencies and thus allow potentially lower drug dosages. However, a proper balancing of compound reactivity is key for the design of covalent binders, to achieve high levels of target inhibition while minimizing promiscuous covalent binding to nontarget proteins. In this work, we demonstrated the possibility to apply the electrophilicity index concept to estimate covalent compound reactivity. We tested this approach on acrylamides, one of the most prominent classes of covalent warheads. Our study clearly demonstrated that, for compounds with molecular weight (MW) below 250 Da, the electrophilicity index can be directly used to estimate compound reactivity. On the other hand, for leadlike molecules (MW > 250 Da) we implemented a new truncation algorithm that has to be applied before reactivity calculations. This algorithm can ensure the localization of HOMO/LUMO orbitals on the compound warhead and thus a correct estimation of its reactivity. Our results also indicate that caution should be used when employing the electrophilicity index to estimate the reactivity of nonterminal acrylamides. The nonparametric nature of this method and its reasonable computational cost make it a suitable tool to support covalent drug design.

Identifying Functional Hotspot Residues for Biased Ligand Design in G-Protein-Coupled Receptors.

G-protein-coupled receptors (GPCRs) mediate multiple signaling pathways in the cell, depending on the agonist that activates the receptor and multiple cellular factors. Agonists that show higher potency to specific signaling pathways over others are known as "biased agonists" and have been shown to have better therapeutic index. Although biased agonists are desirable, their design poses several challenges to date. The number of assays to identify biased agonists seems expensive and tedious. Therefore, computational methods that can reliably calculate the possible bias of various ligands ahead of experiments and provide guidance, will be both cost and time effective. In this work, using the mechanism of allosteric communication from the extracellular region to the intracellular transducer protein coupling region in GPCRs, we have developed a computational method to calculate ligand bias ahead of experiments. We have validated the method for several ß-arrestin-biased agonists in ß2-adrenergic receptor (ß2AR), serotonin receptors 5-HT1B and 5-HT2B and for G-protein-biased agonists in the κ-opioid receptor. Using this computational method, we also performed a blind prediction followed by experimental testing and showed that the agonist carmoterol is ß-arrestin-biased in ß2AR. Additionally, we have identified amino acid residues in the biased agonist binding site in both ß2AR and κ-opioid receptors that are involved in potentiating the ligand bias. We call these residues functional hotspots, and they can be used to derive pharmacophores to design biased agonists in GPCRs.

Multiple Binding Sites Contribute to the Mechanism of Mixed Agonistic and Positive Allosteric Modulators of the Cannabinoid CB1 Receptor.

The cannabinoid CB1 receptor (CB1R) is an abundant metabotropic G-protein-coupled receptor that has been difficult to address therapeutically because of CNS side effects exerted by orthosteric drug candidates. Recent efforts have focused on developing allosteric modulators that target CB1R. Compounds from the recently discovered class of mixed agonistic and positive allosteric modulators (Ago-PAMs) based on 2-phenylindoles have shown promising functional and binding properties as CB1R ligands. Here, we identify binding modes of both the CP 55,940 agonist and GAT228, a 2-phenylindole allosteric modulator, by using our metadynamics simulation protocol, and quantify their affinity and cooperativity by atomistic simulations. We demonstrate the involvement of multiple adjunct binding sites in the Ago-PAM characteristics of the 2-phenylindole modulators and explain their ability to compete with orthosteric agonists at higher concentrations. We validate these results experimentally by showing the contribution of multiple sites on the allosteric binding of ZCZ011, another homologous member of the class, together with the orthosteric agonist.

Development of a 13C NMR Chemical Shift Prediction Procedure Using B3LYP/cc-pVDZ and Empirically Derived Systematic Error Correction Terms: A Computational Small Molecule Structure Elucidation Method.

An accurate and efficient procedure was developed for performing 13C NMR chemical shift calculations employing density functional theory with the gauge invariant atomic orbitals (DFT-GIAO). Benchmarking analysis was carried out, incorporating several density functionals and basis sets commonly used for prediction of 13C NMR chemical shifts, from which the B3LYP/cc-pVDZ level of theory was found to provide accurate results at low computational cost. Statistical analyses from a large data set of 13C NMR chemical shifts in DMSO are presented with TMS as the calculated reference and with empirical scaling parameters obtained from a linear regression analysis. Systematic errors were observed locally for key functional groups and carbon types, and correction factors were determined. The application of this process and associated correction factors enabled assignment of the correct structures of therapeutically relevant compounds in cases where experimental data yielded inconclusive or ambiguous results. Overall, the use of B3LYP/cc-pVDZ with linear scaling and correction terms affords a powerful and efficient tool for structure elucidation.

Energy decomposition analysis approaches and their evaluation on prototypical protein-drug interaction patterns.

The partitioning of the energy in ab initio quantum mechanical calculations into its chemical origins (e.g., electrostatics, exchange-repulsion, polarization, and charge transfer) is a relatively recent development; such concepts of isolating chemically meaningful energy components from the interaction energy have been demonstrated by variational and perturbation based energy decomposition analysis approaches. The variational methods are typically derived from the early energy decomposition analysis of Morokuma [Morokuma, J. Chem. Phys., 1971, 55, 1236], and the perturbation approaches from the popular symmetry-adapted perturbation theory scheme [Jeziorski et al., Methods and Techniques in Computational Chemistry: METECC-94, 1993, ch. 13, p. 79]. Since these early works, many developments have taken place aiming to overcome limitations of the original schemes and provide more chemical significance to the energy components, which are not uniquely defined. In this review, after a brief overview of the origins of these methods we examine the theory behind the currently popular variational and perturbation based methods from the point of view of biochemical applications. We also compare and discuss the chemical relevance of energy components produced by these methods on six test sets that comprise model systems that display interactions typical of biomolecules (such as hydrogen bonding and π-π stacking interactions) including various treatments of the dispersion energy.

SKINK: a web server for string kernel based kink prediction in α-helices.

MOTIVATION: The reasons for distortions from optimal α-helical geometry are widely unknown, but their influences on structural changes of proteins are significant. Hence, their prediction is a crucial problem in structural bioinformatics. Here, we present a new web server, called SKINK, for string kernel based kink prediction. Extending our previous study, we also annotate the most probable kink position in a given α-helix sequence. AVAILABILITY AND IMPLEMENTATION: The SKINK web server is freely accessible at http://biows-inf.zdv.uni-mainz.de/skink. Moreover, SKINK is a module of the BALL software, also freely available at www.ballview.org.

Rodent selectivity of piperidine-4-yl-1H-indoles, a series of CC chemokine receptor-3 (CCR3) antagonists: insights from a receptor model.

Rodent selectivity data of piperidine-4-yl-1H-indoles, a series of CC chemokine receptor-3 (CCR3) antagonists, are presented and discussed as part of an overall optimization effort within this lead compound class. Although attachment of an acidic moiety to the 1-position of the indole led to an overall balanced in vitro profile, in particular reducing inhibition of the hERG channel, potency on the rat and mouse receptor worsened. These findings could be rationalized in the context of a CCR3 homology model.

Density functional theory calculations on entire proteins for free energies of binding: application to a model polar binding site.

In drug optimization calculations, the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method can be used to compute free energies of binding of ligands to proteins. The method involves the evaluation of the energy of configurations in an implicit solvent model. One source of errors is the force field used, which can potentially lead to large errors due to the restrictions in accuracy imposed by its empirical nature. To assess the effect of the force field on the calculation of binding energies, in this article we use large-scale density functional theory (DFT) calculations as an alternative method to evaluate the energies of the configurations in a "QM-PBSA" approach. Our DFT calculations are performed with a near-complete basis set and a minimal parameter implicit solvent model, within the self-consistent calculation, using the ONETEP program on protein-ligand complexes containing more than 2600 atoms. We apply this approach to the T4-lysozyme double mutant L99A/M102Q protein, which is a well-studied model of a polar binding site, using a set of eight small aromatic ligands. We observe that there is very good correlation between the MM and QM binding energies in vacuum but less so in the solvent. The relative binding free energies from DFT are more accurate than the ones from the MM calculations, and give markedly better agreement with experiment for six of the eight ligands. Furthermore, in contrast to MM-PBSA, QM-PBSA is able to correctly predict a nonbinder.

GPCR structures in drug design, emerging opportunities with new structures.

In recent years, GPCR targets from diverse regions of phylogenetic space have been determined. This effort has culminated this year in the determination of representatives of all major classes of GPCRs (A, B, C, and F). Although much of the now well established knowledge on GPCR structures has been known for some years, the new high-resolution structures allow structural insight into the causes of ligand efficacy, biased signaling, and allosteric modulation. In this digest the structural basis for GPCR signaling in the light of the new structures is reviewed and the use of the new non-class A GPCRs for drug design is discussed.

Computational method to identify druggable binding sites that target protein-protein interactions.

Protein-protein interactions are implicated in the pathogenesis of many diseases and are therefore attractive but challenging targets for drug design. One of the challenges in development is the identification of potential druggable binding sites in protein interacting interfaces. Identification of interface surfaces can greatly aid rational drug design of small molecules inhibiting protein-protein interactions. In this work, starting from the structure of a free monomer, we have developed a ligand docking based method, called "FindBindSite" (FBS), to locate protein-protein interacting interface regions and potential druggable sites in this interface. FindBindSite utilizes the results from docking a small and diverse library of small molecules to the entire protein structure. By clustering regions with the highest docked ligand density from FBS, we have shown that these high ligand density regions strongly correlate with the known protein-protein interacting surfaces. We have further predicted potential druggable binding sites on the protein surface using FBS, with druggability being defined as the site with high density of ligands docked. FBS shows a hit rate of 71% with high confidence and 93% with lower confidence for the 41 proteins used for predicting druggable binding sites on the protein-protein interface. Mining the regions of lower ligand density that are contiguous with the high scoring high ligand density regions from FBS, we were able to map 70% of the protein-protein interacting surface in 24 out of 41 structures tested. We also observed that FBS has limited sensitivity to the size and nature of the small molecule library used for docking. The experimentally determined hotspot residues for each protein-protein complex cluster near the best scoring druggable binding sites identified by FBS. These results validate the ability of our technique to identify druggable sites within protein-protein interface regions that have the maximal possibility of interface disruption.

Heteroaromatic π-stacking energy landscapes.

In this study we investigate π-stacking interactions of a variety of aromatic heterocycles with benzene using dispersion corrected density functional theory. We calculate extensive potential energy surfaces for parallel-displaced interaction geometries. We find that dispersion contributes significantly to the interaction energy and is complemented by a varying degree of electrostatic interactions. We identify geometric preferences and minimum interaction energies for a set of 13 5- and 6-membered aromatic heterocycles frequently encountered in small drug-like molecules. We demonstrate that the electrostatic properties of these systems are a key determinant for their orientational preferences. The results of this study can be applied in lead optimization for the improvement of stacking interactions, as it provides detailed energy landscapes for a wide range of coplanar heteroaromatic geometries. These energy landscapes can serve as a guide for ring replacement in structure-based drug design.

Estimation of Electrostatic Interaction Energies on a Trapped-Ion Quantum Computer.

We present the first hardware implementation of electrostatic interaction energies by using a trapped-ion quantum computer. As test system for our computation, we focus on the reduction of NO to N2O catalyzed by a nitric oxide reductase (NOR). The quantum computer is used to generate an approximate ground state within the NOR active space. To efficiently measure the necessary one-particle density matrices, we incorporate fermionic basis rotations into the quantum circuit without extending the circuit length, laying the groundwork for further efficient measurement routines using factorizations. Measurements in the computational basis are then used as inputs for computing the electrostatic interaction energies on a classical computer. Our experimental results strongly agree with classical noise-less simulations of the same circuits, finding electrostatic interaction energies within chemical accuracy despite hardware noise. This work shows that algorithms tailored to specific observables of interest, such as interaction energies, may require significantly fewer quantum resources than individual ground state energies would require in the straightforward supermolecular approach.

Economical and accurate protocol for calculating hydrogen-bond-acceptor strengths.

A series of density functional/basis set combinations and second-order Møller-Plesset calculations have been used to test their ability to reproduce the trends observed experimentally for the strengths of hydrogen-bond acceptors in order to identify computationally efficient techniques for routine use in the computational drug-design process. The effects of functionals, basis sets, counterpoise corrections, and constraints on the optimized geometries were tested and analyzed, and recommendations (M06-2X/cc-pVDZ and X3LYP/cc-pVDZ with single-point counterpoise corrections or X3LYP/aug-cc-pVDZ without counterpoise) were made for suitable moderately high-throughput techniques.

BIreactive: Expanding the Scope of Reactivity Predictions to Propynamides.

We present the first comprehensive study on the prediction of reactivity for propynamides. Covalent inhibitors like propynamides often show improved potency, selectivity, and unique pharmacologic properties compared to their non-covalent counterparts. In order to achieve this, it is essential to tune the reactivity of the warhead. This study shows how three different in silico methods can predict the in vitro properties of propynamides, a covalent warhead class integrated into approved drugs on the market. Whereas the electrophilicity index is only applicable to individual subclasses of substitutions, adduct formation and transition state energies have a good predictability for the in vitro reactivity with glutathione (GSH). In summary, the reported methods are well suited to estimate the reactivity of propynamides. With this knowledge, the fine tuning of the reactivity is possible which leads to a speed up of the design process of covalent drugs.

Towards holistic Compound Quality Scores: Extending ligand efficiency indices with compound pharmacokinetic characteristics.

The suitability of small molecules as oral drugs is often assessed by simple physicochemical rules, the application of ligand efficiency scores or by composite scores based on physicochemical compound properties. These rules and scores are empirical and typically lack mechanistic background, such as information on pharmacokinetics (PK). We introduce new types of Compound Quality Scores (CQS, specifically called dose scores and cmax scores), which explicitly include predicted or, when available, experimental PK parameters and combine these with on-target potency. These CQS scores are surrogates for an estimated dose and corresponding cmax and allow prioritizing of compounds within test cascades as well as before synthesis. We demonstrate the complementarity and, in most cases, superior performance relative to existing efficiency metrics by project examples.