The SAMPL6 challenge on predicting aqueous pKa values from EC-RISM theory

Machine learning methods for pKa prediction of small molecules: Advances and challenges

The acid-base dissociation constant (pK_a) is a fundamental property influencing many ADMET properties of small molecules. However, rapid and accurate pK_a prediction remains a great challenge. In this review, we outline the current advances in machine-learning-based QSAR models for pK_a prediction, including descriptor-based and graph-based approaches, and summarize their pros and cons. Moreover, we highlight the current challenges and future directions regarding experimental data, crucial factors influencing pK_a and in silico prediction tools. We hope that this review can provide a practical guidance for the follow-up studies.

Improving Small Molecule pKa Prediction Using Transfer Learning With Graph Neural Networks

Enumerating protonation states and calculating microstate pK_a values of small molecules is an important yet challenging task for lead optimization and molecular modeling. Commercial and non-commercial solutions have notable limitations such as restrictive and expensive licenses, high CPU/GPU hour requirements, or the need for expert knowledge to set up and use. We present a graph neural network model that is trained on 714,906 calculated microstate pK_a predictions from molecules obtained from the ChEMBL database. The model is fine-tuned on a set of 5,994 experimental pK_a values significantly improving its performance on two challenging test sets. Combining the graph neural network model with Dimorphite-DL, an open-source program for enumerating ionization states, we have developed the open-source Python package pkasolver, which is able to generate and enumerate protonation states and calculate pK_a values with high accuracy.

Implementation and Optimization of the Embedded Cluster Reference Interaction Site Model with Atomic Charges

In this work, we implemented the embedded cluster reference interaction site model (EC-RISM) originally developed by Kloss, Heil, and Kast (J. Phys. Chem. B2008, 112, 4337–4343). This method combines quantum mechanical calculations with the 3D reference interaction site model (3D-RISM). Numerous options, such as buffer, grid space, basis set, charge model, water model, closure relation, and so forth, were investigated to find the best settings. Additionally, the small point charges, which are derived from the solvent distribution from the 3D-RISM solution to represent the solvent in the QM calculation, were neglected to reduce the overhead without the loss of accuracy. On the MNSOL[a], MNSOL, and FreeSolv databases, our implemented and optimized method provides solvation free energies in water with 5.70, 6.32, and 6.44 kJ/mol root-mean-square deviations, respectively, but with different settings, 5.22, 6.08, and 6.63 kJ/mol can also be achieved. Only solvent models containing fitting parameters, like COSMO-RS and EC-RISM with universal correction and directly used electrostatic potential, perform better than our EC-RISM implementation with atomic charges.

Recent Developments of Computational Methods for pKa Prediction Based on Electronic Structure Theory with Solvation Models

The protonation/deprotonation reaction is one of the most fundamental processes in solutions and biological systems. Compounds with dissociative functional groups change their charge states by protonation/deprotonation. This change not only significantly alters the physical properties of a compound itself, but also has a profound effect on the surrounding molecules. In this paper, we review our recent developments of the methods for predicting the Ka, the equilibrium constant for protonation reactions or acid dissociation reactions. The pKa, which is a logarithm of Ka, is proportional to the reaction Gibbs energy of the protonation reaction, and the reaction free energy can be determined by electronic structure calculations with solvation models. The charge of the compound changes before and after protonation; therefore, the solvent effect plays an important role in determining the reaction Gibbs energy. Here, we review two solvation models: the continuum model, and the integral equation theory of molecular liquids. Furthermore, the reaction Gibbs energy calculations for the protonation reactions require special attention to the handling of dissociated protons. An efficient method for handling the free energy of dissociated protons will also be reviewed.

Asymmetric Interplay Between K+ and Blocker and Atomistic Parameters From Physiological Experiments Quantify K+ Channel Blocker Release

Modulating the activity of ion channels by blockers yields information on both the mode of drug action and on the biophysics of ion transport. Here we investigate the interplay between ions in the selectivity filter (SF) of K⁺ channels and the release kinetics of the blocker tetrapropylammonium in the model channel Kcv_NTS. A quantitative expression calculates blocker release rate constants directly from voltage-dependent ion occupation probabilities in the SF. The latter are obtained by a kinetic model of single-channel currents recorded in the absence of the blocker. The resulting model contains only two adjustable parameters of ion-blocker interaction and holds for both symmetric and asymmetric ionic conditions. This data-derived model is corroborated by 3D reference interaction site model (3D RISM) calculations on several model systems, which show that the K⁺ occupation probability is unaffected by the blocker, a direct consequence of the strength of the ion-carbonyl attraction in the SF, independent of the specific protein background. Hence, Kcv_NTS channel blocker release kinetics can be reduced to a small number of system-specific parameters. The pore-independent asymmetric interplay between K⁺ and blocker ions potentially allows for generalizing these results to similar potassium channels.

A Joint Venture of Ab Initio Molecular Dynamics, Coupled Cluster Electronic Structure Methods, and Liquid-State Theory to Compute Accurate Isotropic Hyperfine Constants of Nitroxide Probes in Water

The isotropic hyperfine coupling constant (HFCC, Aiso) of a pH-sensitive spin probe in a solution, HMI (2,2,3,4,5,5-hexamethylimidazolidin-1-oxyl, C9H19N2O) in water, is computed using an ensemble of state-of-the-art computational techniques and is gauged against X-band continuous wave electron paramagnetic resonance (EPR) measurement spectra at room temperature. Fundamentally, the investigation aims to delineate the cutting edge of current first-principles-based calculations of EPR parameters in aqueous solutions based on using rigorous statistical mechanics combined with correlated electronic structure techniques. In particular, the impact of solvation is described by exploiting fully atomistic, RISM integral equation, and implicit solvation approaches as offered by ab initio molecular dynamics (AIMD) of the periodic bulk solution (using the spin-polarized revPBE0-D3 hybrid functional), embedded cluster reference interaction site model integral equation theory (EC-RISM), and polarizable continuum embedding (using CPCM) of microsolvated complexes, respectively. HFCCs are obtained from efficient coupled cluster calculations (using open-shell DLPNO-CCSD theory) as well as from hybrid density functional theory (using revPBE0-D3). Re-solvation of "vertically desolvated" spin probe configuration snapshots by EC-RISM embedding is shown to provide significantly improved results compared to CPCM since only the former captures the inherent structural heterogeneity of the solvent close to the spin probe. The average values of the Aiso parameter obtained based on configurational statistics using explicit water within AIMD and from EC-RISM solvation are found to be satisfactorily close. Using either such explicit or RISM solvation in conjunction with DLPNO-CCSD calculations of the HFCCs provides an average Aiso parameter for HMI in aqueous solution at 300 K and 1 bar that is in good agreement with the experimentally determined one. The developed computational strategy is general in the sense that it can be readily applied to other spin probes of similar molecular complexity, to aqueous solutions beyond ambient conditions, as well as to other solvents in the longer run.

SAMPL7 physical property prediction from EC-RISM theory

Inspired by the successful application of the embedded cluster reference interaction site model (EC-RISM), a combination of quantum–mechanical calculations with three-dimensional RISM theory to predict Gibbs energies of species in solution within the SAMPL6.1 (acidity constants, pK_a) and SAMPL6.2 (octanol–water partition coefficients, log P) the methodology was applied to the recent SAMPL7 physical property challenge on aqueous pK_a and octanol–water log P values. Not part of the challenge but provided by the organizers, we also computed distribution coefficients log D_7.4 from predicted pK_a and log P data. While macroscopic pK_a predictions compared very favorably with experimental data (root mean square error, RMSE 0.72 pK units), the performance of the log P model (RMSE 1.84) fell behind expectations from the SAMPL6.2 challenge, leading to reasonable log D_7.4 predictions (RMSE 1.69) from combining the independent calculations. In the post-submission phase, conformations generated by different methodology yielded results that did not significantly improve the original predictions. While overall satisfactory compared to previous log D challenges, the predicted data suggest that further effort is needed for optimizing the robustness of the partition coefficient model within EC-RISM calculations and for shaping the agreement between experimental conditions and the corresponding model description.

Evaluation of log P, pKa, and log D predictions from the SAMPL7 blind challenge

The Statistical Assessment of Modeling of Proteins and Ligands (SAMPL) challenges focuses the computational modeling community on areas in need of improvement for rational drug design. The SAMPL7 physical property challenge dealt with prediction of octanol-water partition coefficients and pK_a for 22 compounds. The dataset was composed of a series of N-acylsulfonamides and related bioisosteres. 17 research groups participated in the log P challenge, submitting 33 blind submissions total. For the pK_a challenge, 7 different groups participated, submitting 9 blind submissions in total. Overall, the accuracy of octanol-water log P predictions in the SAMPL7 challenge was lower than octanol-water log P predictions in SAMPL6, likely due to a more diverse dataset. Compared to the SAMPL6 pK_a challenge, accuracy remains unchanged in SAMPL7. Interestingly, here, though macroscopic pK_a values were often predicted with reasonable accuracy, there was dramatically more disagreement among participants as to which microscopic transitions produced these values (with methods often disagreeing even as to the sign of the free energy change associated with certain transitions), indicating far more work needs to be done on pK_a prediction methods.

Quantum-mechanical property prediction of solvated drug molecules: what have we learned from a decade of SAMPL blind prediction challenges?

Joint academic-industrial projects supporting drug discovery are frequently pursued to deploy and benchmark cutting-edge methodical developments from academia in a real-world industrial environment at different scales. The dimensionality of tasks ranges from small molecule physicochemical property assessment over protein-ligand interaction up to statistical analyses of biological data. This way, method development and usability both benefit from insights gained at both ends, when predictiveness and readiness of novel approaches are confirmed, but the pharmaceutical drug makers get early access to novel tools for the quality of drug products and benefit of patients. Quantum-mechanical and simulation methods particularly fall into this group of methods, as they require skills and expense in their development but also significant resources in their application, thus are comparatively slowly dripping into the realm of industrial use. Nevertheless, these physics-based methods are becoming more and more useful. Starting with a general overview of these and in particular quantum-mechanical methods for drug discovery we review a decade-long and ongoing collaboration between Sanofi and the Kast group focused on the application of the embedded cluster reference interaction site model (EC-RISM), a solvation model for quantum chemistry, to study small molecule chemistry in the context of joint participation in several SAMPL (Statistical Assessment of Modeling of Proteins and Ligands) blind prediction challenges. Starting with early application to tautomer equilibria in water (SAMPL2) the methodology was further developed to allow for challenge contributions related to predictions of distribution coefficients (SAMPL5) and acidity constants (SAMPL6) over the years. Particular emphasis is put on a frequently overlooked aspect of measuring the quality of models, namely the retrospective analysis of earlier datasets and predictions in light of more recent and advanced developments. We therefore demonstrate the performance of the current methodical state of the art as developed and optimized for the SAMPL6 pKa and octanol-water log P challenges when re-applied to the earlier SAMPL5 cyclohexane-water log D and SAMPL2 tautomer equilibria datasets. Systematic improvement is not consistently found throughout despite the similarity of the problem class, i.e. protonation reactions and phase distribution. Hence, it is possible to learn about hidden bias in model assessment, as results derived from more elaborate methods do not necessarily improve quantitative agreement. This indicates the role of chance or coincidence for model development on the one hand which allows for the identification of systematic error and opportunities toward improvement and reveals possible sources of experimental uncertainty on the other. These insights are particularly useful for further academia-industry collaborations, as both partners are then enabled to optimize both the computational and experimental settings for data generation.

Overview of the SAMPL6 pKa challenge: evaluating small molecule microscopic and macroscopic pKa predictions

The prediction of acid dissociation constants (pKa) is a prerequisite for predicting many other properties of a small molecule, such as its protein-ligand binding affinity, distribution coefficient (log D), membrane permeability, and solubility. The prediction of each of these properties requires knowledge of the relevant protonation states and solution free energy penalties of each state. The SAMPL6 pKa Challenge was the first time that a separate challenge was conducted for evaluating pKa predictions as part of the Statistical Assessment of Modeling of Proteins and Ligands (SAMPL) exercises. This challenge was motivated by significant inaccuracies observed in prior physical property prediction challenges, such as the SAMPL5 log D Challenge, caused by protonation state and pKa prediction issues. The goal of the pKa challenge was to assess the performance of contemporary pKa prediction methods for drug-like molecules. The challenge set was composed of 24 small molecules that resembled fragments of kinase inhibitors, a number of which were multiprotic. Eleven research groups contributed blind predictions for a total of 37 pKa distinct prediction methods. In addition to blinded submissions, four widely used pKa prediction methods were included in the analysis as reference methods. Collecting both microscopic and macroscopic pKa predictions allowed in-depth evaluation of pKa prediction performance. This article highlights deficiencies of typical pKa prediction evaluation approaches when the distinction between microscopic and macroscopic pKas is ignored; in particular, we suggest more stringent evaluation criteria for microscopic and macroscopic pKa predictions guided by the available experimental data. Top-performing submissions for macroscopic pKa predictions achieved RMSE of 0.7-1.0 pKa units and included both quantum chemical and empirical approaches, where the total number of extra or missing macroscopic pKas predicted by these submissions were fewer than 8 for 24 molecules. A large number of submissions had RMSE spanning 1-3 pKa units. Molecules with sulfur-containing heterocycles or iodo and bromo groups were less accurately predicted on average considering all methods evaluated. For a subset of molecules, we utilized experimentally-determined microstates based on NMR to evaluate the dominant tautomer predictions for each macroscopic state. Prediction of dominant tautomers was a major source of error for microscopic pKa predictions, especially errors in charged tautomers. The degree of inaccuracy in pKa predictions observed in this challenge is detrimental to the protein-ligand binding affinity predictions due to errors in dominant protonation state predictions and the calculation of free energy corrections for multiple protonation states. Underestimation of ligand pKa by 1 unit can lead to errors in binding free energy errors up to 1.2 kcal/mol. The SAMPL6 pKa Challenge demonstrated the need for improving pKa prediction methods for drug-like molecules, especially for challenging moieties and multiprotic molecules.

Predicting pKa Using a Combination of Semi-Empirical Quantum Mechanics and Radial Basis Function Methods

The acid dissociation constant (pK_a) has an important influence on molecular properties crucial to compound development in synthesis, formulation, and optimization of absorption, distribution, metabolism, and excretion properties. We will present a method that combines quantum mechanical calculations, at a semi-empirical level of theory, with machine learning to accurately predict pK_a for a diverse range of mono- and polyprotic compounds. The resulting model has been tested on two external data sets, one specifically used to test pK_a prediction methods (SAMPL6) and the second covering known drugs containing basic functionalities. Both sets were predicted with excellent accuracy (root-mean-square errors of 0.7-1.0 log units), comparable to other methodologies using a much higher level of theory and computational cost.

SAMPL6 logP challenge: machine learning and quantum mechanical approaches

Two different types of approaches: (a) approaches that combine quantitative structure activity relationships, quantum mechanical electronic structure methods, and machine-learning and, (b) electronic structure vertical solvation approaches, were used to predict the logP coefficients of 11 molecules as part of the SAMPL6 logP blind prediction challenge. Using electronic structures optimized with density functional theory (DFT), several molecular descriptors were calculated for each molecule, including van der Waals areas and volumes, HOMO/LUMO energies, dipole moments, polarizabilities, and electrophilic and nucleophilic superdelocalizabilities. A multilinear regression model and a partial least squares model were used to train a set of 97 molecules. As well, descriptors were generated using the molecular operating environment and used to create additional machine learning models. Electronic structure vertical solvation approaches considered include DFT and the domain-based local pair natural orbital methods combined with the solvated variant of the correlation consistent composite approach.

Assessing the accuracy of octanol-water partition coefficient predictions in the SAMPL6 Part II log P Challenge

The SAMPL Challenges aim to focus the biomolecular and physical modeling community on issues that limit the accuracy of predictive modeling of protein-ligand binding for rational drug design. In the SAMPL5 log D Challenge, designed to benchmark the accuracy of methods for predicting drug-like small molecule transfer free energies from aqueous to nonpolar phases, participants found it difficult to make accurate predictions due to the complexity of protonation state issues. In the SAMPL6 log P Challenge, we asked participants to make blind predictions of the octanol-water partition coefficients of neutral species of 11 compounds and assessed how well these methods performed absent the complication of protonation state effects. This challenge builds on the SAMPL6 p[Formula: see text] Challenge, which asked participants to predict p[Formula: see text] values of a superset of the compounds considered in this log P challenge. Blind prediction sets of 91 prediction methods were collected from 27 research groups, spanning a variety of quantum mechanics (QM) or molecular mechanics (MM)-based physical methods, knowledge-based empirical methods, and mixed approaches. There was a 50% increase in the number of participating groups and a 20% increase in the number of submissions compared to the SAMPL5 log D Challenge. Overall, the accuracy of octanol-water log P predictions in SAMPL6 Challenge was higher than cyclohexane-water log D predictions in SAMPL5, likely because modeling only the neutral species was necessary for log P and several categories of method benefited from the vast amounts of experimental octanol-water log P data. There were many highly accurate methods: 10 diverse methods achieved RMSE less than 0.5 log P units. These included QM-based methods, empirical methods, and mixed methods with physical modeling supported with empirical corrections. A comparison of physical modeling methods showed that QM-based methods outperformed MM-based methods. The average RMSE of the most accurate five MM-based, QM-based, empirical, and mixed approach methods based on RMSE were 0.92 ± 0.13, 0.48 ± 0.06, 0.47 ± 0.05, and 0.50 ± 0.06, respectively.

Tautomeric Equilibria of Nucleobases in the Hachimoji Expanded Genetic Alphabet

Evolution has yielded biopolymers that are constructed from exactly four building blocks and are able to support Darwinian evolution. Synthetic biology aims to extend this alphabet, and we recently showed that 8-letter (hachimoji) DNA can support rule-based information encoding. One source of replicative error in non-natural DNA-like systems, however, is the occurrence of alternative tautomeric forms, which pair differently. Unfortunately, little is known about how structural modifications impact free-energy differences between tautomers of the non-natural nucleobases used in the hachimoji expanded genetic alphabet. Determining experimental tautomer ratios is technically difficult, and so, strategies for improving hachimoji DNA replication efficiency will benefit from accurate computational predictions of equilibrium tautomeric ratios. We now report that high-level quantum-chemical calculations in aqueous solution by the embedded cluster reference interaction site model, benchmarked against free-energy molecular simulations for solvation thermodynamics, provide useful quantitative information on the tautomer ratios of both Watson-Crick and hachimoji nucleobases. In agreement with previous computational studies, all four Watson-Crick nucleobases adopt essentially only one tautomer in water. This is not the case, however, for non-natural nucleobases and their analogues. For example, although the enols of isoguanine and a series of related purines are not populated in water, these heterocycles possess N₁-H and N₃-H keto tautomers that are similar in energy, thereby adversely impacting accurate nucleobase pairing. These robust computational strategies offer a firm basis for improving experimental measurements of tautomeric ratios, which are currently limited to studying molecules that exist only as two tautomers in solution.

Standard state free energies, not pKas, are ideal for describing small molecule protonation and tautomeric states

The pK_a is the standard measure used to describe the aqueous proton affinity of a compound, indicating the proton concentration (pH) at which two protonation states (e.g. A^- and AH) have equal free energy. However, compounds can have additional protonation states (e.g. AH₂⁺), and may assume multiple tautomeric forms, with the protons in different positions (microstates). Macroscopic pK_as give the pH where the molecule changes its total number of protons, while microscopic pK_as identify the tautomeric states involved. As tautomers have the same number of protons, the free energy difference between them and their relative probability is pH independent so there is no pK_a connecting them. The question arises: What is the best way to describe protonation equilibria of a complex molecule in any pH range? Knowing the number of protons and the relative free energy of all microstates at a single pH, ∆G°, provides all the information needed to determine the free energy, and thus the probability of each microstate at each pH. Microstate probabilities as a function of pH generate titration curves that highlight the low energy, observable microstates, which can then be compared with experiment. A network description connecting microstates as nodes makes it straightforward to test thermodynamic consistency of microstate free energies. The utility of this analysis is illustrated by a description of one molecule from the SAMPL6 Blind pK_a Prediction Challenge. Analysis of microstate ∆G°s also makes a more compact way to archive and compare the pH dependent behavior of compounds with multiple protonatable sites.

The SAMPL6 challenge on predicting octanol-water partition coefficients from EC-RISM theory

Results are reported for octanol–water partition coefficients (log P) of the neutral states of drug-like molecules provided during the SAMPL6 (Statistical Assessment of Modeling of Proteins and Ligands) blind prediction challenge from applying the “embedded cluster reference interaction site model” (EC-RISM) as a solvation model for quantum-chemical calculations. Following the strategy outlined during earlier SAMPL challenges we first train 1- and 2-parameter water-free (“dry”) and water-saturated (“wet”) models for n-octanol solvation Gibbs energies with respect to experimental values from the “Minnesota Solvation Database” (MNSOL), yielding a root mean square error (RMSE) of 1.5 kcal mol⁻¹ for the best-performing 2-parameter wet model, while the optimal water model developed for the pK_a part of the SAMPL6 challenge is kept unchanged (RMSE 1.6 kcal mol⁻¹ for neutral compounds from a model trained on both neutral and ionic species). Applying these models to the blind prediction set yields a log P RMSE of less than 0.5 for our best model (2-parameters, wet). Further analysis of our results reveals that a single compound is responsible for most of the error, SM15, without which the RMSE drops to 0.2. Since this is the only compound in the challenge dataset with a hydroxyl group we investigate other alcohols for which Gibbs energy of solvation data for both water and n-octanol are available in the MNSOL database to demonstrate a systematic cause of error and to discuss strategies for improvement.

High pressure response of 1H NMR chemical shifts of purine nucleotides

The study of the pressure response by NMR spectroscopy provides information on the thermodynamics of conformational equilibria in proteins and nucleic acids. For obtaining a database for expected pressure effects on free nucleotides and nucleotides bound in macromolecular complexes, the pressure response of ¹H chemical shifts and J-coupling constants of the purine 5'-ribonucleotides AMP, ADP, ATP, GMP, GDP, and GTP were studied in the absence and presence of Mg²⁺-ions. Experiments are supported by quantum-chemical calculations of populations and chemical shift differences in order to corroborate structural interpretations and to estimate missing data for AMP. The preference of the ribose S puckering obtained from the analysis of the experimental J-couplings is also confirmed by the calculations. In addition, the pressure response of the non-hydrolysable GTP analogues GppNHp, GppCH₂p, and GTPγS was examined within a pressure range up to 200 MPa. As observed earlier for ³¹P NMR chemical shifts of these nucleotides the pressure dependence of chemical shifts is clearly non-linear in most cases. In di- and tri-phospho nucleosides, the resonances of the two protons bound to the ribose 5' carbon are non-equivalent and can be observed separately. The gg-rotamer at C4'- C5' bond is strongly preferred and the downfield shifted resonance can be assigned to the H5″ proton in the nucleotides. In contrast, in adenosine itself the frequencies of the two resonances are interchanged.

pKa calculations for tautomerizable and conformationally flexible molecules: partition function vs. state transition approach

Calculations of acidities of molecules with multiple tautomeric and/or conformational states require adequate treatment of the relative energetics of accessible states accompanied by a statistical-mechanical formulation of their contribution to the macroscopic pK_a value. Here, we demonstrate rigorously the formal equivalence of two such approaches: a partition function treatment and statistics over transitions between molecular tautomeric and conformational states in the limit of a theory that does not require adjustment by empirical parameters correcting energetic values. However, for a frequently employed correction scheme, linear scaling of (free) energies and regression with respect to reference data taking an additive constant into account, this equivalence breaks down if more than one acid or base state is involved. The consequences of the resulting inconsistency are discussed on our datasets developed for aqueous pK_a predictions during the recent SAMPL6 challenge, where molecular state energetics were computed based on the "embedded cluster reference interaction site model" (EC-RISM). This method couples integral equation theory as a solvation model to quantum-chemical calculations and yielded a test set root mean square error of 1.1 pK units from a partition function ansatz. For all practical purposes, the present results indicate that a state transition approach yields comparable accuracy despite the formal theoretical inconsistency, and that an additive regression intercept, which is strictly constant in the limit of large compound mass only, is a valid approximation. Graphical abstract Embedded cluster reference interaction site model-derived vs. experimental pK_a for the test set calculated with either the partition function (blue) or the state transition approach (red), using m as a free parameter.