SAMPL6: Calculation of macroscopic pK_a values from ab initio quantum mechanical free energies

1. Introduction

The SAMPL (Statistical Assessment of the Modeling of Proteins and Ligands) challenges allow the molecular modeling community to assess, in “blind” conditions, the accuracy and efficiency of current computational chemistry methods and tools, leading to continuous improvements of the available computational methods. The previous SAMPL challenges [1–5] involved hydration free energy calculations, with the exception of the last edition, SAMPL5, which was dedicated to the prediction of distribution coefficients [6]. Our past participations in SAMPL challenges [7–9] represented unique opportunities for us to test our approaches and to develop and improve new computational tools. In 2018, the SAMPL6 challenge focused on the prediction of microscopic and macroscopic pK_a values for fragment-like organic compounds.

The equilibrium acid dissociation reaction in aqueous solution

with acid dissociation constant $K_{\text{a}}\text{=}\left[{\text{A}}^{-}\right]\left[{\text{H}}^{+}\right]\text{/}\left[\text{HA}\right]$ is of broad importance in biological systems, in synthetic chemistry, and pharmacology [10–14]. The pK_a, defined as

for the standard state concentration c₀ = 1 mol/l, measures thermodynamic acidity. The theoretical prediction of pK_a values is an ongoing challenge [15]. In the narrow realm of protein biochemistry, good progress has been made in calculating the physiologically important changes in pK_as of standard amino acid residues in different environments with accuracies better than 1 pK_a unit [12], especially with constant pH molecular dynamics simulations [16–19], which have been applied to study a wide range of phenomena [20–22]. Absolute pK_a calculations of arbitrary molecules using physics-based quantum chemistry approaches (as opposed to machine learning (ML) ones) have been more challenging and accuracy of 1 pK_a unit has been difficult to achieve consistently [15, 23] whereas a range of methods can achieve “chemical accuracy” (defined as 2.5 pK_a units by Ho and Coote [15]). The clear advantage of ab initio approaches is that they can be applied to any novel compound. Here we report on pK_a calculations of the 24 compounds in the SAMPL6 challenge (Fig. 1) with a quantum-chemical approach originally developed by Muckerman et al [24]. The SAMPL6 compounds are, however, chemically more complex and 23 contain multiple titratable protons and, in some cases, tautomers so that macroscopic pK_a have to be calculated.

The calculation of microscopic pK_as, i.e., the free energy difference for the deprotonation reaction Eq. 1 at the standard state (concentration 1 mol/l and temperature T = 298.15 K, indicated by the superscript “*”)

is straightforward using quantum chemical gas-phase calculations. However, it is well-known [15, 23] that direct calculations lead to large errors in the calculated pK_as, mainly due to the poor continuum solvation models that have to be employed in order to obtain free energies in solution. One approach to correct for these systematic errors is to generate a model to correct the raw quantum chemical free energies [24]. We generated linear models from a training set with 38 simple compounds with experimentally known pK_a (Fig. 2 and 3). We fit a global model to all the data (the global linear model) and we split the training set with a simple classifier, namely the charge of the acid, yielding a piecewise linear model with separate linear functions for neutral and cationic acids. We calculated the macroscopic pK_as for all 24 SAMPL6 compounds and compared the accuracy of the three approaches [QM computed (raw), linear fit global, and linear fit piecewise].

2. Methods

Following Muckerman et al [24], our strategy was to compute gas-phase free energy differences

(denoted as standard state free energies at 1 atm pressure and 298.15 K) for the deprotonation reaction for all titratable protons,

To obtain solution free energy differences corresponding to Eq. 1,

(where the standard state refers to 1 mol/l), a solvation free energy contribution $\Delta G_{\text{solv}}^{\text{ο}}$ is added to the gas-phase free energies of the acid HA and the base A⁻ from Eq. 4,

with $\Delta G^{\circ \to *}=1.894\text{ kcal/mol}$ accounting for the change in standard state in the gas phase. The free energy of the proton in the gas phase is calculated analytically in the ideal gas limit (the Sackur-Tetrode equation [25]), $G^{\text{ο}}\left(H_{\text{(g)}}^{+}\right)\text{=}-\text{6}.28\text{kcal/mol},$ and for the solvation free energy of the proton we chose the same value as Muckerman et al [24], G*(HA_(s)) = −272.2 kcal/mol although other values are also discussed in the literature [15, 26]. With $\Delta G_{\text{(s)}}^{\text{*}}$ , the pK_a is calculated from Eq. 3.

As described in detail in Section 2.2, the directly calculated pK_a values have fairly poor accuracy and thus we derive a simple linear estimator to correct for shortcomings in the solvation model [24]. The linear model is based on our own training data set (described in the next section) and the resulting estimator ℒ is applied to the pK_a from Eq. 3 to obtain improved predictions for the SAMPL6 data set, $\text{p}{K}_{\text{a}}\text{=}ℒ\text{[p}{K}_{\text{a}}^{\text{calc}}\text{]}$.

2.1. Data sets

The QM1 subset of the training set contains 21 neutral acids belonging to several chemical families (Fig. 2): mono- (1), di- (2) and tri- (3) protic inorganic acids, aliphatic (4) and aromatic (5) sulfonic acids, diversely substituted carboxylic acids (6–11) and alcohols (12–17), phenols (18 and 19), phthalimide (20) and uracil (21). The experimental pK_a values of these compounds range from −3.00 to 17.10 (Table 1). The QM2 subset contains 17 compounds that are cationic acids (Fig. 3): hydrazine (22), guanidine (23), aliphatic mono- (24), di- (25) and tri- (26) substituted amines, diversely substituted aromatic amines (27–31) and pyridines (32–38). These compounds possess experimental pK_a values from 0.49 to 13.60 (Table 1).

The SAMPL6 data set consisted of 24 fragment-like small organic molecules (Fig. 1) with unknown pK_a values that were selected for their similarity to kinase inhibitors and for experimental tractability. It was provided by the SAMPL6 organizers through the public repository https://github.com/MobleyLab/SAMPL6 as computer-generated microstates in SMILES format. The protonation state for each microstate was computed with an in-house script using the C ACT VS Chemoinfor- matics Toolkit [27] (Xemistry GmbH, https://www.xemistry.com/), allowing the classification of microstates in two groups, neutral acids and cationic acids, for which different correction factors were applied in the approach using the piecewise linear model.

Three-dimensional coordinates for all compounds were generated in MOL2 format using CORINA version 3.60 (http://www.molecular-networks.com), then converted into the Gaussian input format using an in-house script. The PDF 3d files, which can be visualized with Adobe Acrobat Reader (https://get.adobe.com/fr/reader/) were generated with Cactvs.

2.2. Quantum chemical microscopic pK_a calculations

Gas-phase geometry optimization and frequency calculation of the protonated and deprotonated forms were performed at the B3LYP/6–311+G(d,p) level using Gaussian 09 version D.01 [28] to obtain $\Delta G_{\text{(s)}}^{\text{ο}}$. A single-point free energy evaluation at the same level using the Conductor-like Polarizable Continuum Solvation Model (CPCM) [29–32] and UAHF radii as implemented in Gaussian 09 version D.01 [28] yielded the solvation free energy $\Delta G_{\text{solv}}^{\text{ο}}$ so that $\Delta G_{\text{(s)}}^{\text{*}}$ (Eq. 7) and an estimate for the pK_a associated with this protonation/deprotonation event could be calculated via Eqs. 6 and 3.

In some cases, the geometry optimization did not converge with Gaussian 09 version D.01, but was successful with the version A.02 of Gaussian 09. Geometry optimization for microstates SM04_micro016, SM07_micro016, SM17_micro008 and SM17_micro009 did not converge in any conditions.

Muckerman et al [24] recognized systematic errors related to the solvation contribution $\Delta G_{\text{solv}}^{\text{*}}$ as responsible for poor accuracy, namely the solvation model under-solvates weak acids and over-solvates strong acids. They proposed a physically-motivated correction

to $\Delta G_{\text{solv}}^{°}$ with the linear model

The parameters a₀ and a₁ are determined from a training set by linear regression. In order to apply the correction Eq. 9 to compounds with unknown pK_a, a linear estimator ℒ can be derived by substituting $\text{p}{K}_{\text{a}}^{\exp }\approx \text{p}{K}_{\text{a}}^{\text{calc}}\text{+ }\Delta G_{\text{corr}}^{\text{*}}/(RT\ln \text{10})$ in Eq. 9 and solving for $\Delta G_{\text{corr}}^{\text{*}}$ to yield

The linear estimator ℒ with parameters a₀ and a₁ for the microscopic pK_a is

2.3. Microstates vs Macrostates

We consider each tautomer of the acid HA and the base A⁻ as a microstate with label i. The set of microstates with the same total number of protons N_i = N is labeled the macrostate N. The macroscopic pK_a characterizes the transitions between any of the microstates with N protons to any microstate with N − 1 protons.

In general, the free energy difference between two states (micro or macro states) that are separated by a single protonation process (i.e., the free energy to go from N to N − 1 associated protons) is

where P(N − 1) and P(N) are the probabilities of observing the system with N − 1 and N associated protons respectively and β = (RT)⁻¹. The Henderson-Hasselbalch equation

can be rewritten in terms of the free energy of protonation $\Delta G_{N-1,N}$ (Eq. 12) to give

2.4. Calculation of macroscopic pK_as from microscopic pK_as

The microscopic pK_a values correspond to free energy differences $\Delta G_{ij}(\text{pH})\text{ = }{G}_j(\text{pH})-G_i\text{(pH)}$ between microstates i and j (Eq. 14b); for notational convenience we drop the explicit pH dependence in the following for all free energies. Each state has a pH-dependent associated free energy G_i, which is not known. Constructing the G_i from the differences between them is not straightforward because these calculated free energy differences come with unknown errors that prevent, for example, that the sum along any closed thermodynamic cycle $i\to j\to k\to \cdot \cdot \cdot \to i$ is exactly zero as required by the fact that the G_i are thermodynamic state functions. We construct a set of M microstate free energies ${\left\{G_i\right\}}_{i=1}^M$ that is most consistent with the calculated (“measured”) {ΔG_ij} using a maximum-likelihood estimator [33] based on the likelihood function

where we assumed normal distribution of errors with constant standard deviation. The product runs over all pairs (i, j) for which calculated ΔGij are available. L is proportional to the probability Ρ({ΔG_ij}|{G_i}) that we could observe the measured data (all the calculated Δ G_ij) if we were given a specific set of the G_i (our model parameters). Maximizing the log-likelihood ln L (using functions in SciPy [34]) as a function of all the G_i provides the set ${\left\{G_i\right\}}_{i=1}^M$ that is most consistent with the given measurements {Δ Gij}. Further details and more general applications of this approach will be published elsewhere (I.M. Kenney et al, in preparation).

In order to calculate the macroscopic pK_as, we begin by calculating the free energy of protonation using principles of equilibrium statistical mechanics [25]. The probability of observing a macrostate with N associated protons is

where the sum is over all accessible microstates with free energy G_i, δ_Νi,_N is unity when the microstate i has N protons and null otherwise, and Z is the partition function, defined by

Eq. 16 combined with the general expression for the free energy of protonation (Eq. 12) yields the effective macroscopic protonation free energy as a function of the G_i,

ΔG_N-1,N is a function of the pH of the system and the microscopic pK_as relevant to the macrostate N. Together with Eq. 14a, Eq. 18 allows us to calculate the macroscopic pK_a value for removing the N^th proton from a molecule. With all microstate free energies ${\left\{G_i\right\}}_{i=1}^M$ known for a given pH value it is also straightforward to compute the pH-dependent microstate probabilities

where all terms depend on pH.

2.5. Error analysis

The difference between experimental and computed pK_a values (“signed error”) for each compound, labeled with its identification code ‘id’, was calculated as

The root-mean-square error (RMSE) was determined from the individual errors Δ as

the mean absolute error (MAE) as

and the signed mean error (ME, also called the “mean signed error”, MSE) as

We also report the Pearson correlation coefficient R² and the slope m of a linear regression to the data, as computed with the function scipy.stats.linregress( ) in the SciPy package [34].

The quantum chemical single point free energy calculations do not have a statistical error and we have not yet implemented the calculation of an error bound in the maximum likelihood estimator for the G_i. Therefore, all pK_a are provided without a statistical error. Judging from the performance of the training data set and the post- hoc analysis of the SAMPL6 compounds (see Results), the accuracy of the calculated pK_a values is 1–2 pK_a units.

Calculated pK_a were compared to experimental values with the script typeIII_analysis.py as provided by the SAMPL6 organizers in the public repository https://github.com/MobleyLab/SAMPL6. Calculated values were matched to experimental ones with the Hungarian algorithm, which finds the optimum pairing between two sets by minimizing the linear sum of squared errors.

3. Results and Discussion

4. Conclusions

Compared to other methods in the SAMPL6 challenge, our approach has below- average accuracy (Fig. 9 and Figs. S4–S7 in the Electronic Supplementary Material) and its computational cost is also higher than ML-based approaches (not considering the cost for compiling and validating the data and training the ML model). A key advantage of our approach is its generality as it does not depend on training on specific data sets although below we note that the quality of the training set for the correction step is a possible concern. With the linear model, which was derived from a very small and simple training set (38 compounds), we remove some of the errors related to the QM method used and its implementation in Gaussian (e.g., the implicit solvation model). The quality of the prediction is mostly independent of the structure, i.e., it can predict organic compounds from different families and even inorganic compounds with similar level of accuracy. In comparison, purely ML-based methods are trained on large experimental data sets (containing several thousands or tens of thousands compounds) and they can be vulnerable to chemical families that are not represented in the training set. Our approach appears reasonably robust because for our training set we obtain the same slope on the global data set and on the individual subsets, which are chemically quite different. The results for the SAMPL6 data set are also fairly insensitive to the fitting approach used (global or piecewise linear model), further demonstrating robustness. The correlations with experimental data are generally good but suffer from systematic errors, possibly from differences between the training set and the SAMPL6 set that bias the estimator that is needed to correct the raw QM pK_a values. The statistical measures indicate clear room from improvement. It appears that a better correction scheme, using a larger data set that better matches the test data set with respect to its hydrogen bonding properties and is generally more representative of drug-like molecules could improve the predictions, perhaps in conjunction with more sophisticated classifiers and estimators than simple separation by charge and linear regression. We expect that improvements in the model physics, namely in the treatment of solvation, could also lead to further increases in accuracy.

We currently consider the method described here (and originally developed by Muckerman et al [24]) as an acceptable compromise between speed, accuracy and generality across the chemical space. It seems especially useful when one encounters novel compounds and wants to assess them based on their absolute pK_a values. The calculations are tractable with typical computational resources, absolute pKas are accurate to about 2 units (within the “chemical accuracy” range [15]) and do not seem to be biased with respect to specific chemical groups, and thus the relative ordering of compounds is also meaningful.

Supplementary Material