General Theory of Fragment Linking in Molecular Design: Why Fragment Linking Rarely Succeeds and How to Improve Outcomes.

Linking two fragments binding in nearby subpockets together has become an important technique in fragment-based drug discovery to optimize the binding potency of fragment hits. Despite the expected favorable translational and orientational entropic contribution to the binding free energy of the linked molecule, brute force enumeration of chemical linker for linking fragments is rarely successful, and the vast majority of linked molecules do not exhibit the expected gains of binding potency. In this paper, we examine the physical factors that contribute to the change of binding free energy from fragment linking and develop a method to rigorously calculate these different physical contributions. We find from these analyses that multiple confounding factors make successful fragment linking strategies rare, including (1) possible change of the binding mode of the fragments in the linked state compared to separate binding of the fragments, (2) unfavorable intramolecular strain energy of the bioactive conformation of the linked molecule, (3) unfavorable interaction between the linker and the protein, (4) favorable interaction energies between two fragments in solution when not chemically linked that offset the expected entropy loss for the formation of fragment pair, (5) complex compensating configurational entropic effects beyond the simplistic rotational and translational analysis. We here have applied a statistically mechanically rigorous approach to compute the fragment linking coefficients of 10 pharmaceutically interesting systems and quantify the contribution of each physical component to the binding free energy of the linked molecule. Based on these studies, we have found that the change in the relative configurational entropy of the two fragments in the protein binding pocket (a term neglected to our knowledge in all previous analyses) substantially offsets the favorable expected rotational and translational entropic contributions to the binding free energy of the linked molecule. This configurational restriction of the fragments in the binding pocket of the proteins is found to be, in our analysis, the dominant reason why most fragment linking strategies do not exhibit the expected gains of binding potency. These findings have further provided rich physical insights, which we expect should facilitate more successful fragment linking strategies to be formulated in the future.

Pattern-free generation and quantum mechanical scoring of ring-chain tautomers.

In contrast to the computational generation of conventional tautomers, the analogous operation that would produce ring-chain tautomers is rarely available in cheminformatics codes. This is partly due to the perceived unimportance of ring-chain tautomerism and partly because specialized algorithms are required to realize the non-local proton transfers that occur during ring-chain rearrangement. Nevertheless, for some types of organic compounds, including sugars, warfarin analogs, fluorescein dyes and some drug-like compounds, ring-chain tautomerism cannot be ignored. In this work, a novel ring-chain tautomer generation algorithm is presented. It differs from previously proposed solutions in that it does not rely on hard-coded patterns of proton migrations and bond rearrangements, and should therefore be more general and maintainable. We deploy this algorithm as part of a workflow which provides an automated solution for tautomer generation and scoring. The workflow identifies protonatable and deprotonatable sites in the molecule using a previously described approach based on rapid micro-pKa prediction. These data are used to distribute the active protons among the protonatable sites exhaustively, at which point alternate resonance structures are considered to obtain pairs of atoms with opposite formal charge. These pairs are connected with a single bond and a 3D undistorted geometry is generated. The scoring of the generated tautomers is performed with a subsequent density functional theory calculation employing an implicit solvent model. We demonstrate the performance of our workflow on several types of organic molecules known to exist in ring-chain tautomeric equilibria in solution. In particular, we show that some ring-chain tautomers not found using previously published algorithms are successfully located by ours.

Advancing Free-Energy Calculations of Metalloenzymes in Drug Discovery via Implementation of LFMM Potentials.

To address some of the inherent challenges in modeling metalloenzymes, we here report an extension to the functional form of the OPLS3e force field to include terms adopted from the ligand field molecular mechanics (LFMM) model, including the angular overlap and Morse potential terms. The integration of these terms with OPLS3e, herein referred to as OPLS3e+M, improves the description of metal-ligand interactions and provides accurate relative binding energies and geometric preferences of transition-metal complexes by training to gas-phase density functional theory (DFT) energies. For [Cu(H2O)4]2+, OPLS3e+M significantly improves H2O binding energies and the geometric preference of the tetra-aqua Cu2+ complex. In addition, we conduct free-energy perturbation calculations on two pharmaceutically relevant metalloenzyme targets, which include chemical modifications at varying proximity to the binding-site metals, including changes to the metal-binding moiety of the ligand itself. The extensions made to OPLS3e lead to accurate predicted relative binding free energies for these series (mean unsigned error of 1.29 kcal mol-1). Our results provide evidence that integration of the LFMM model with OPLS3e can be utilized to predict thermodynamic quantities for such systems near chemical accuracy. With these improvements, we anticipate that robust free-energy perturbation calculations can be employed to accelerate the drug development efforts for metalloenzyme targets.

Evaluation of Free Energy Calculations for the Prioritization of Macrocycle Synthesis.

A tremendous research and development effort was exerted toward combating chronic hepatitis C, ultimately leading to curative oral treatments, all of which are targeting viral proteins. Despite the advantage of numerous targets allowing for broad hepatitis C virus (HCV) genotype coverage, the only host target inhibitors that advanced into clinical development were Cyclosporin A based cyclophilin inhibitors. While cyclosporin-based molecules typically require a fermentation process, Gilead successfully pursued a fully synthetic, oral program based on Sanglifehrin A. The drug discovery process, though greatly helped by facile crystallography, was still hampered by the limitations in the accuracy of predictive computational methods for prioritizing compound ideas. Recent advances in accuracy and speed of free energy perturbation (FEP) methods, however, are attractive for prioritizing and derisking synthetically challenging molecules and potentially could have had a significant impact on the speed of the development of this program. Here in our simulated prospective study, the binding free energies of 26 macrocyclic cyclophilin inhibitors were blindly predicted using FEP+ to test this hypothesis. The predictions had a low mean unsigned error (MUE) (1.1 kcal/mol) and accurately reproduced many design decisions from the program, suggesting that FEP+ has the potential to drive synthetic chemistry efforts by more accurately ranking compounds with nonintuitive structure-activity relationships (SARs).

Toward Atomistic Modeling of Irreversible Covalent Inhibitor Binding Kinetics.

Covalent inhibitors have emerged as an important drug class in recent years, largely due to their many unique advantages as compared to noncovalent inhibitors, including longer duration of action, lower prolonged systemic exposure, higher potency, and selectivity. However, the potential off-target toxicity of covalent inhibitors, particularly of irreversible covalent inhibitors, represents a great challenge in covalent drug development. Therefore, accurate calculation of protein covalent inhibitor reaction kinetics to guide the design of selective inhibitors would greatly benefit covalent drug discovery efforts. In the present paper, we present a computational method to calculate the relative reaction kinetics between congeneric irreversible covalent inhibitors and their protein receptors. The method combines density functional theory calculations of the transition state barrier height of the rate-limiting step for reaction between the warhead of the inhibitor and a single protein residue, and molecular-mechanics-based free energy calculations to account for the interactions between the ligand in the transition state and the protein environment. The method was tested on four pharmaceutically interesting irreversible covalent binding systems involving 28 ligands; the mean unsigned error (MUE) of the relative reaction rate for all pairs of ligands between the predictions and experimental results for these tested systems is 0.79 log unit. This is to our knowledge the first time where the reaction kinetics of protein irreversible covalent inhibition have been directly calculated with physics-based free energy calculation methods and transition state theory. We anticipate the outstanding accuracy demonstrated here across a broad range of target classes will have a strong impact on the design of selective covalent inhibitors.

Generation of Tautomers Using Micro-p Ka's.

Solutions of organic molecules containing one or more heterocycles with conjugated bonds may exist as a mixture of tautomers, but typically only a few of them are significantly populated even though the potential number grows combinatorially with the number of protonation and deprotonation sites. Generating the most stable tautomers from a given input structure is an important and challenging task, and numerous algorithms to tackle it have been proposed in the literature. This work describes a novel approach for tautomer prediction that involves the combined use of molecular mechanics, semiempirical quantum chemistry, and density functional theory. The key idea in our method is to identify the protonation and deprotonation sites using estimated micro-p Ka's for every atom in the molecule as well as in its nearest protonated and deprotonated forms. To generate tautomers in a systematic way with minimal bias, we then consider the full set of tautomers that arise from the combinatorial distribution of all such mobile protons among all protonatable sites, with efficient postprocessing to screen away high-energy species. To estimate the micro-p Ka's, we present a new method designed for the current task, but we emphasize that any alternative method can be used in conjunction with our basic algorithm. Our approach is therefore grounded in the computational prediction of physical properties in aqueous solution, in contrast to other approaches that may rely on the use of hard-coded rules of proton distribution, previously observed tautomerization patterns from a known chemical space, or human input. We present examples of the application of our algorithm to organic and drug-like molecules, with a focus on novel structures where traditional methods are expected to perform worse.

Rigorous Free Energy Perturbation Approach to Estimating Relative Binding Affinities between Ligands with Multiple Protonation and Tautomeric States.

Accurate prediction of ligand binding affinities is of key importance in small molecule lead optimization and a central task in computational medicinal chemistry. Over the years, advances in both computer hardware and computational methodologies have established free energy perturbation (FEP) methods as among the most reliable and rigorous approaches to compute protein-ligand binding free energies. However, accurate description of ionization and tautomerism of ligands is still a major challenge in structure-based prediction of binding affinities. Druglike molecules are often weak acid or bases with multiple accessible protonation and tautomeric states that can contribute significantly to the binding process. To address this issue, we introduce in this work the p Ka and tautomeric state correction approach. This approach is based on free energy perturbation formalism and provides a rigorous treatment of the ionization and tautomeric equilibria of ligands in solution and in the protein complexes. A series of Kinesin Spindle Protein (KSP) and Factor Xa inhibitor molecules were used as test cases. Our results demonstrate that the p Ka and tautomeric state correction approach is able to rigorously and accurately incorporate multiple protonation and tautomeric states in the binding affinity calculations.

On the Upper Limits of Oxidation States in Chemistry.

The concept of oxidation state (OS) is based on the concept of Lewis electron pairs, in which the bonding electrons are assigned to the more electronegative element. This approach is useful for keeping track of the electrons, predicting chemical trends, and guiding syntheses. Experimental and quantum-chemical results reveal a limit near +8 for the highest OS in stable neutral chemical substances under ambient conditions. OS=+9 was observed for the isolated [IrO4 ]+ cation in vacuum. The prediction of OS=+10 for isolated [PtO4 ]2+ cations is confirmed computationally for low temperatures only, but hasn't yet been experimentally verified. For high OS species, oxidation of the ligands, for example, of O-2 with formation of . O-1 and O-O bonds, and partial reduction of the metal center may be favorable, possibly leading to non-Lewis type structures.

Weighted Averaging Scheme and Local Atomic Descriptor for pKa Prediction Based on Density Functional Theory.

As a continuation of our work on developing a density functional theory-based pKa predictor, we present conceptual improvements to our previously published shell model, which is a hierarchical organization of pKa training sets and which, in principle, covers all chemical space. The improvements concern the way the studied chemical compound is associated with the data points from the training sets. By introducing a new descriptor of the local atomic environment which foregoes dependence on chemical bonding and connectivity, we are able to automatically locate molecules from the training set that are most relevant to the proton dissociation equilibrium under study. This new scheme leads to the prediction of a single pKa value weighted across multiple training sets and thus patches a defect disclosed in the formulation of our previous model. Using the new parametrization approach, the pKa prediction gets rid of outliers reported in previous applications of our approach, eliminates ambiguity in interpreting the results, and improves the overall accuracy. Our new treatment accounts for multiple conformations both on the level of energetics and parametrization. Illustrative results are shown for several types of chemical structures containing guanidine, amidine, amine, and phenol functional groups, and which are representative of practically important large and flexible drug-like molecules. Our method's performance is compared to the performance of other previously published pKa prediction methods. Further possible improvements to the organization of the training sets and the potential application of our new local atomic descriptor to other kinds of parametrizations are discussed.

Accurate and Reliable Prediction of the Binding Affinities of Macrocycles to Their Protein Targets.

Macrocycles have been emerging as a very important drug class in the past few decades largely due to their expanded chemical diversity benefiting from advances in synthetic methods. Macrocyclization has been recognized as an effective way to restrict the conformational space of acyclic small molecule inhibitors with the hope of improving potency, selectivity, and metabolic stability. Because of their relatively larger size as compared to typical small molecule drugs and the complexity of the structures, efficient sampling of the accessible macrocycle conformational space and accurate prediction of their binding affinities to their target protein receptors poses a great challenge of central importance in computational macrocycle drug design. In this article, we present a novel method for relative binding free energy calculations between macrocycles with different ring sizes and between the macrocycles and their corresponding acyclic counterparts. We have applied the method to seven pharmaceutically interesting data sets taken from recent drug discovery projects including 33 macrocyclic ligands covering a diverse chemical space. The predicted binding free energies are in good agreement with experimental data with an overall root-mean-square error (RMSE) of 0.94 kcal/mol. This is to our knowledge the first time where the free energy of the macrocyclization of linear molecules has been directly calculated with rigorous physics-based free energy calculation methods, and we anticipate the outstanding accuracy demonstrated here across a broad range of target classes may have significant implications for macrocycle drug discovery.

Revised M06-L functional for improved accuracy on chemical reaction barrier heights, noncovalent interactions, and solid-state physics.

We present the revM06-L functional, which we designed by optimizing against a larger database than had been used for Minnesota 2006 local functional (M06-L) and by using smoothness restraints. The optimization strategy reduced the number of parameters from 34 to 31 because we removed some large terms that increased the required size of the quadrature grid and the number of self-consistent-field iterations. The mean unsigned error (MUE) of revM06-L on 422 chemical energies is 3.07 kcal/mol, which is improved from 3.57 kcal/mol calculated by M06-L. The MUE of revM06-L for the chemical reaction barrier height database (BH76) is 1.98 kcal/mol, which is improved by more than a factor of 2 with respect to the M06-L functional. The revM06-L functional gives the best result among local functionals tested for the noncovalent interaction database (NC51), with an MUE of only 0.36 kcal/mol, and the MUE of revM06-L for the solid-state lattice constant database (LC17) is half that for M06-L. The revM06-L functional also yields smoother potential curves, and it predicts more-accurate results than M06-L for seven out of eight diversified test sets not used for parameterization. We conclude that the revM06-L functional is well suited for a broad range of applications in chemistry and condensed-matter physics.

Perspective: Kohn-Sham density functional theory descending a staircase.

This article presents a perspective on Kohn-Sham density functional theory (KS-DFT) for electronic structure calculations in chemical physics. This theory is in widespread use for applications to both molecules and solids. We pay special attention to several aspects where there are both concerns and progress toward solutions. These include: 1. The treatment of open-shell and inherently multiconfigurational systems (the latter are often called multireference systems and are variously classified as having strong correlation, near-degeneracy correlation, or high static correlation; KS-DFT must treat these systems with broken-symmetry determinants). 2. The treatment of noncovalent interactions. 3. The choice between developing new functionals by parametrization, by theoretical constraints, or by a combination. 4. The ingredients of the exchange-correlation functionals used by KS-DFT, including spin densities, the magnitudes of their gradients, spin-specific kinetic energy densities, nonlocal exchange (Hartree-Fock exchange), nonlocal correlation, and subshell-dependent corrections (DFT+U). 5. The quest for a universal functional, where we summarize some of the success of the latest Minnesota functionals, namely MN15-L and MN15, which were obtained by optimization against diverse databases. 6. Time-dependent density functional theory, which is an extension of DFT to treat time-dependent problems and excited states. The review is a snapshot of a rapidly moving field, and-like Marcel Duchamp-we hope to convey progress in a stimulating way.

Oxidation State 10 Exists.

In a recent paper, Wang et al. found an iridium-containing compound with a formal oxidation state of 9. This is the highest oxidation state ever found in a stable compound. To learn if this is the highest chemical oxidation state possible, Kohn-Sham density functional theory was used to study various compounds, including PdO4 (2+) , PtO4 (2+) , PtO3 F2 (2+) , PtO4 OH(+) , PtO5 , and PtO4 SH(+) , in which the metal has an oxidation state of 10. It was found that PtO4 (2+) has a metastable state that is kinetically stable with a barrier height for decomposition of 31 kcal mol(-1) and a calculated lifetime of 0.9 years. All other compounds studied would readily decompose to lower oxidation states.

MN15-L: A New Local Exchange-Correlation Functional for Kohn-Sham Density Functional Theory with Broad Accuracy for Atoms, Molecules, and Solids.

Kohn-Sham density functional theory is widely used for applications of electronic structure theory in chemistry, materials science, and condensed-matter physics, but the accuracy depends on the quality of the exchange-correlation functional. Here, we present a new local exchange-correlation functional called MN15-L that predicts accurate results for a broad range of molecular and solid-state properties including main-group bond energies, transition metal bond energies, reaction barrier heights, noncovalent interactions, atomic excitation energies, ionization potentials, electron affinities, total atomic energies, hydrocarbon thermochemistry, and lattice constants of solids. The MN15-L functional has the same mathematical form as a previous meta-nonseparable gradient approximation exchange-correlation functional, MN12-L, but it is improved because we optimized it against a larger database, designated 2015A, and included smoothness restraints; the optimization has a much better representation of transition metals. The mean unsigned error on 422 chemical energies is 2.32 kcal/mol, which is the best among all tested functionals, with or without nonlocal exchange. The MN15-L functional also provides good results for test sets that are outside the training set. A key issue is that the functional is local (no nonlocal exchange or nonlocal correlation), which makes it relatively economical for treating large and complex systems and solids. Another key advantage is that medium-range correlation energy is built in so that one does not need to add damped dispersion by molecular mechanics in order to predict accurate noncovalent binding energies. We believe that the MN15-L functional should be useful for a wide variety of applications in chemistry, physics, materials science, and molecular biology.

Correction: MN15: A Kohn-Sham global-hybrid exchange-correlation density functional with broad accuracy for multi-reference and single-reference systems and noncovalent interactions.

[This corrects the article DOI: 10.1039/C6SC00705H.].

MN15: A Kohn-Sham global-hybrid exchange-correlation density functional with broad accuracy for multi-reference and single-reference systems and noncovalent interactions.

Kohn-Sham density functionals are widely used; however, no currently available exchange-correlation functional can predict all chemical properties with chemical accuracy. Here we report a new functional, called MN15, that has broader accuracy than any previously available one. The properties considered in the parameterization include bond energies, atomization energies, ionization potentials, electron affinities, proton affinities, reaction barrier heights, noncovalent interactions, hydrocarbon thermochemistry, isomerization energies, electronic excitation energies, absolute atomic energies, and molecular structures. When compared with 82 other density functionals that have been defined in the literature, MN15 gives the second smallest mean unsigned error (MUE) for 54 data on inherently multiconfigurational systems, the smallest MUE for 313 single-reference chemical data, and the smallest MUE on 87 noncovalent data, with MUEs for these three categories of 4.75, 1.85, and 0.25 kcal mol-1, respectively, as compared to the average MUEs of the other 82 functionals of 14.0, 4.63, and 1.98 kcal mol-1. The MUE for 17 absolute atomic energies is 7.4 kcal mol-1 as compared to an average MUE of the other 82 functionals of 34.6 kcal mol-1. We further tested MN15 for 10 transition-metal coordination energies, the entire S66x8 database of noncovalent interactions, 21 transition-metal reaction barrier heights, 69 electronic excitation energies of organic molecules, 31 semiconductor band gaps, seven transition-metal dimer bond lengths, and 193 bond lengths of 47 organic molecules. The MN15 functional not only performs very well for our training set, which has 481 pieces of data, but also performs very well for our test set, which has 823 data that are not in our training set. The test set includes both ground-state properties and molecular excitation energies. For the latter MN15 achieves simultaneous accuracy for both valence and Rydberg electronic excitations when used with linear-response time-dependent density functional theory, with an MUE of less than 0.3 eV for both types of excitations.

Nonseparable exchange-correlation functional for molecules, including homogeneous catalysis involving transition metals.

The goal of this work is to develop a gradient approximation to the exchange-correlation functional of Kohn-Sham density functional theory for treating molecular problems with a special emphasis on the prediction of quantities important for homogeneous catalysis and other molecular energetics. Our training and validation of exchange-correlation functionals is organized in terms of databases and subdatabases. The key properties required for homogeneous catalysis are main group bond energies (database MGBE137), transition metal bond energies (database TMBE32), reaction barrier heights (database BH76), and molecular structures (database MS10). We also consider 26 other databases, most of which are subdatabases of a newly extended broad database called Database 2015, which is presented in the present article and in its ESI. Based on the mathematical form of a nonseparable gradient approximation (NGA), as first employed in the N12 functional, we design a new functional by using Database 2015 and by adding smoothness constraints to the optimization of the functional. The resulting functional is called the gradient approximation for molecules, or GAM. The GAM functional gives better results for MGBE137, TMBE32, and BH76 than any available generalized gradient approximation (GGA) or than N12. The GAM functional also gives reasonable results for MS10 with an MUE of 0.018 Å. The GAM functional provides good results both within the training sets and outside the training sets. The convergence tests and the smooth curves of exchange-correlation enhancement factor as a function of the reduced density gradient show that the GAM functional is a smooth functional that should not lead to extra expense or instability in optimizations. NGAs, like GGAs, have the advantage over meta-GGAs and hybrid GGAs of respectively smaller grid-size requirements for integrations and lower costs for extended systems. These computational advantages combined with the relatively high accuracy for all the key properties needed for molecular catalysis make the GAM functional very promising for future applications.