A Weighted Averaging Scheme and a Local Atomic Descriptor for pKa

Subscriber access provided by READING UNIV

Article

A Weighted Averaging Scheme and a Local Atomic Descriptor for pKa Prediction Based on Density Functional Theory Haoyu S, Yu, Mark A. Watson, and Art D. Bochevarov J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.7b00537 • Publication Date (Web): 22 Jan 2018 Downloaded from http://pubs.acs.org on January 24, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Information and Modeling is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Information and Modeling

A Weighted Averaging Scheme and a Local Atomic Descriptor for pKa Prediction Based on Density Functional Theory Haoyu S. Yu, Mark A. Watson, and Art D. Bochevarov Schrödinger, Inc., 120 West 45th St, New York, New York 10036

Abstract 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 parameterization 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 parameterization. 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

1

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 50

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 parameterizations are discussed.

1

INTRODUCTION

For small and medium-sized organic molecules, the computation of the negative logarithm of the proton dissociation constant, pKa , based on the density functional theory (DFT) 1–5 can be viewed as complimentary to the other popular approach, an empirical one. 6–8 Both approaches find their own niches of application. The former describes protonation processes as well as conformational and steric effects explicitly and, due to its moderate dependence on parameterization, promises higher accuracy for novel types of structures. However, the high computational cost of DFT-based methods does not permit their application to more than a handful of structures in a single project. For a discussion of challenges facing quantum chemical pKa prediction of complex and conformationally flexible organic molecules see Ref. 9 The empirical approach’s best asset is its high computational efficiency yielding an almost instant pKa prediction. This advantage can be used to screen large libraries of compounds, with the accuracy deemed to be high for the structures that are similar to those from the underlying training set, and less reliable for unusual structures not well represented in it. The empirical approach, directly exploiting structure-property correlations, might also provide an insufficient insight into the physical and chemical phenomena responsible for the predicted pKa values. Some pKa transitions observed in a titration experiment is a result of the protonation or deprotonation of multiple functional groups. This can happen when the molecule has several locations involved in protonation or deprotonation equilibria around the same pH. It may also be the case that the chemical substance studied in a titration experiment undergoes tautomerization or even chemically reacts with the solvent. The single, formal chemical formula ascribed to the substance will thus not be representative of the chemical processes actually responsible for its measured pKa . An empirical pKa prediction method which associates the chemical formula with

2


Page 3 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


its observed pKa might not reveal the complexity of the dissociation processes actually occurring in solution, even if the empirical prediction method turns out to be accurate. In a previous publication 5 we described an automated workflow which uses a combined molecular mechanics (MM) and DFT approach to predict pKa for a large part of the practically important chemical space of small and medium-sized molecules. The general organization of the previously described workflow, as well as its modified version presented in this work, are very similar and can be summarized as follows. The workflow acts on a 3D structure in which one of the atoms is marked by the user, which is necessary to define the specific protonation or deprotonation process to model. Following a combination of gas phase and solution phase calculations on the protonated and the deprotonated structures, with an option to generate conformers thereof, the algorithm uses the computed energetics data to assemble the free energy of the modeled process. The free energy is further adjusted for the symmetry of the molecule and to account for specific interactions with the solvent, which produces a final (micro-) pKa value that, in case of e single protonation site, should be comparable to the experimentally measured pKa . See Section 5 for more detail regarding the choice of conformations, the application of parameterization, and combining multiple micro-pKa ’s into macro-pKa where necessary. The MM part of the previously published algorithm was used for taking care of conformational flexibility, which is absolutely necessary for accurate pKa predictions of large organic molecules. Our workflow differs from other DFT approaches reported in the literature in that it works on an extensive area of practically important chemical space, includes a sophisticated selection and treatment of conformers, and can act on large and flexible molecules producing results of high accuracy. In addition, the implementation is automated to a considerable degree, requiring only a single structure and a specification of the protonation/deprotonation process as input. One deficiency of our method was the dependence of the generated pKa values on the choice of the underlying parameterization (or training) test set, such that at the end of a calculation typically several pKa values would be reported, each associated with its own training set. Because we automatically select training sets on the basis of the notion of chemical similarity, 10 more

3



than one training set, and therefore more than one parameterization model, can be regarded as appropriate for describing the pKa of a given compound. For instance, a structure of the type RN = C(NHR0 )(NHR00 ), assumed to be protonatable at its sp2 -hybridized nitrogen atom, might be viewed as a guanidine, amidine, partly substituted amidine, or even imine, depending on how many covalent bonds in the immediate vicinity of the protonated atom one considers. Consequently, for the above-mentioned structure, one might compute several pKa values, all corresponding to different parameterization models based on training sets of guanidines, amidines, etc. A human researcher usually has little difficulty associating a given chemical compound with a broader class of compounds, typically characterized by their functional group which predominantly contributes to the chemical phenomenon of interest, like a protonation equilibrium. The human designation of the functional group is, however, subjective, and while it does not appear questionable for clear-cut cases, it becomes tantalizingly difficult for borderline chemical structures. Take phenols for an example. Popular educational resources define phenols as compounds in which the hydroxy-group is connected to either the benzene ring 11 or more generally to an aromatic ring or arene, 12–15 although in the latter case all sources discuss the benzene ring as the only representative of the aromatic rings. The definition through aromatic rings formally hinges on the vaguely defined concept of aromaticity 16 and suggests that naphtols, hydroxyquinolines, hydroxypyridines, and even the unstable 2-hydroxyfuran could be formally recognized as phenols. Most chemists would probably draw the line between phenols and other classes of hydroxy-compounds somewhere in the middle of Figure 1, thinking of the extreme case of 2-hydroxyfuran as an enol, due to its well-known tautomeric equilibrium strongly shifted to the keto-form. 17 This essentially semantic issue unfortunately carries over to other areas of research and engenders problems in the specification of parameters in models that are supposed to describe a broad range of chemical compounds. A hydroxyl-substituted, presumably aromatic heterocycle the tautomeric equilibria of which have not been well studied would present a challenging case indeed: can parameters developed for hydroxy-benzenes (phenols proper) be used for this heterocycle or not? Our experience of developing parameters for pKa prediction indicates that the subjective association of a given compound or its

4


Page 4 of 50

Page 5 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


fragment with a broader class of compounds produces artifacts and practical inconveniences: (i) final results depend on the chosen association with a class of compounds due to a different use of parameters; (ii) costly human involvement is required to classify the compound.

Figure 1: The identity of phenols is difficult to pinpoint but is important for assigning model parameters to chemical structures. The image shows five different molecules which can all be thought of as “phenols”, according to the definition that phenols are structures with the hydroxy group attached to an aromatic ring. 12–15 Many chemists would probably draw a (dashed) line between phenols (closer to the left) and non-phenols (closer to the right), as shown in the figure. Thus, it is important to develop an automated parameterization scheme that would alleviate or eliminate the problems mentioned in the previous paragraph. In the following section we introduce a novel descriptor of local chemical environment and the corresponding local similarity measure. The approach based on this descriptor brings us closer to an automatic classification of chemical structures and assigning parameters associated with the structure’s atoms. The descriptor chosen in our work does not depend explicitly on the notion of chemical bonding. As such, it has a number of advantages in comparison with the approaches that rely on connectivity and bonding information. Although we apply the new descriptor to a pKa prediction protocol, we believe it can be used in other computational models that deal with assignment of local parameters. We discuss this possibility in the text. We combine our new similarity measure with a specially constructed distance metric to associate the closeness of a given structure to training sets or classes of chemical compounds. This leads to a formulation of a practical scheme to produce a single-value, weighted pKa prediction and promises several improvements in the way we construct our training sets. These enhancements allow us to eliminate an important problem inherent in the formulation of our previous model, which has been discussed earlier. 5 We apply the new methodology to the calculation of pKa of 5



compounds with two types of functional groups: those that can be viewed as belonging to one class of compounds only (for example, amines), and those that might be regarded as belonging to several classes (for example, guanidines can be regarded as a special case of amidines). These applications demonstrate that the improved model proposed in this paper is capable of eliminating outliers associated with the conceptual problem from the earlier model, and has little effect on the results that were already accurate. We conclude with a discussion of how the proposed local atomic descriptors can be used for a streamlined, on-the-fly construction of pKa parameters, tailored to a given compound and avoiding costly human involvement in classification of chemical compounds and construction of training sets. It might be also possible to use the proposed local atomic descriptor for improving the organization of parameterization in other models where atomic parameters have to be set (e.g. force fields, solvation models). An implementation of these ideas is planned for future work.

2 LOCAL ATOMIC DESCRIPTORS A variety of descriptors representing local atomic environments have been developed and discussed in the literature. 18–21 Such descriptors have been actively used in machine learning methods. 21–23 Below we are going to describe a novel 3D local atomic descriptor suitable for a pKa prediction method that utilizes empirical and DFT data. Because proton dissociation processes are local in nature, for our purposes we seek a descriptor that characterizes the immediate environment of the dissociating proton. We also prefer a descriptor that abandons dependence on assigned chemical bonds between atoms. This is for three reasons. First, there are practical difficulties of comparing structures on the basis of chemical bonding. If chemical bonding is at the foundation of the comparison method, then the user or developer usually has to come up with SMARTS patterns, which are textual, 1D-representations of 2D chemical formulas that explicitly specify bonding between atoms. 24 These representations have proved to be very useful in various cheminformatics applications 25,26 and parameterization schemes. 27,28

6


Page 6 of 50

Page 7 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


However, for a human, SMARTS patterns are not trivial to write and read for anything but the simplest cases. These patterns become particularly complicated if they are meant to target certain structures and exclude other structures; for example, a pattern that describes “phenols excluding ortho-nitro- and ortho-nitrosophenols” (see this pattern in Table S23). Often SMARTS patterns require careful testing, as unintended mistakes in them, easy to make in practice, can result in mismatches of structures and large errors in final results. Entering and maintaining hundreds or even thousands of SMARTS patterns in extensive parameterization models becomes prohibitively expensive. The second reason is that bonding is only well defined for lighter elements and simpler organic structures. In structures of inorganic and organometallic origin the assignment of covalent chemical bonds is particularly ambiguous, with the consequence that models, parameters, or definitions relying on chemical bonds and atom types become fragile or inapplicable. Other languages for textual representations of molecular structures exist, 29,30 and although it would be interesting to evaluate them for our purposes, we suspect that these alternative representations will suffer from similar defects. Finally, connectivity patterns suffer from ambiguities due to multiple possible resonance structures (mesomers). It is important to note that results produced by quantum chemical methods such as DFT normally do not depend on the resonance structure (connectivity, bond orders, and formal atomic charges) entered as input. However, if the quantum chemical method is supplemented with a structure matching empirical component, as is the case of our pKa prediction method, and if this empirical component is sensitive to the choice of the resonance structure, then the overall result might depend on the resonance structure too. There is usually a preferred and universally accepted resonance structure for each functional group, and the parameterization specific to this functional group can be associated with the more common resonance structure. However, it is disconcerting that the result of the calculation and the quality of pKa prediction might depend on the resonance structure chosen by the user. There are of course computational canonicalization schemes which return a “standard” resonance structure for any given one, but they may lead to resonance structures that are at variance with the user’s

7



Page 8 of 50

intuition. The inconsistency between the machine representation and the user preference for the resonance structure may be regarded as a simple inconvenience. But it needs to be verified that models which use resonance structures to assign parameters on atoms and functional groups still produce good quality results if fed machine-generated, canonical resonance structures, particularly in cases when the parameters were trained on simpler resonance structures fixed by the user. Our local atomic descriptor L can be computed for every atom A in the molecule. The descriptor L computed at atom A (or at local environment A in general) will be denoted as LA , whereas the descriptor L computed for molecule X (at an unspecified atom) or even for a set of compounds S will be denoted as L(X) or L(S), respectively. The descriptor consists of two parts. The first part incorporates information about the geometric surroundings of A (the location of other atoms in the vicinity). The second part characterizes the electrostatic effects near A. We can write the local atomic descriptor for atom A as h i ~ A , QA LA = M

(1)

In practice, the information contained in LA should be sufficient to distinguish one local environment from another. ~ , is the multidimensional coordination number. It characterizes The first part of the descriptor, M the geometric environment of the atom and must be “myopic”: it must be “blind” to distant geometric effects, in accordance with the notion of locality of proton dissociation effects. For example, distant conformational effects should not influence the local energetics of a protonation or deprotonation at a certain atom. In constructing this part of the descriptor, we borrow the mathematical ideas of Grimme and co-workers from a different area of research, 31 adjusting these ideas to our purpose. We try, however, to stick to the notation of the original work. The multidimensional coordination number is defined as a vector ~ A = {CNA (H), CNA (He), CNA (Li), ...} M

8


(2)

Page 9 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where each component is computed for an element of the periodic table. In principle, the dimension of M is the number of elements in the periodic table, but in practice it can be truncated to the number of elements in actual use. The definition of CNA (E) is as follows: N (E)

CNA (E) =

X (B∈E)6=A

1 1 + exp [−k1 (k2 (RA + RB )/rAB − 1)]

(3)

The summation goes over all atoms B of element E, which are not atom A, and N (E) is the number of atoms of type E. Just like in the work of Grimme and co-workers, k1 and k2 are adjustable parameters, RA is a scaled covalent single bond radius of atom A, and rAB is the interatomic distance between atoms A and B. The physical meaning of the coordination number can be further elucidated if we examine its mathematical properties. Consider for simplicity formula (3) in which the summation goes over a single atom B 6= A. In the work of Grimme and co-workers, as well as in our work (vide infra) k1 is typically around 16 and k2 is a little larger than 1. Under these conditions, near the equilibrium atomic distance rAB the expression exp [−k1 (k2 (RA + RB )/rAB − 1)] in the denominator will be very small and CNA (E) will be approaching 1. The meaning of this value is that A is connected to one atom B. At small interatomic distances (much smaller than equilibrium rAB ) the contribution converges to 1. At an infinite separation of atoms A and B the coordination number CNA (E) tends to 0, meaning that the atoms are not bonded. At intermediate distances the connectivity of A to B is between 1 and 0. So, when we regain the summation over all atoms B of type (E), the meaning of CNA (E) is the number of atoms of type (E) effectively coordinated (connected) with atom A. For simple cases in which covalent bonding is clearly defined the coordination number closely resembles the sum over the number of covalent connections for corresponding atom type. For example, CNC (H) for an equilibrium structure of methane is close to 4.0, and CNC (C) in every carbon atom in benzene is close to 2.0. One difference between the earlier and the current definition of the coordination number (CN) is that Grimme and co-workers compute CNA without distinguishing the elements with which atom

9



A is coordinated, whereas we split the coordination number into components each characterizing the coordination of atom A with atoms of a particular type. In such a way, our descriptor would be able to tell between CH4 and CCl4 , whereas the descriptor that does not differentiate between the elements would not. So, MA is a generalization of CNA from Ref., 31 and the latter can be obtained from the former by summing up all the components in equation (2). Splitting coordination number into element components allows us to get a finer characterization of the local environment. In some formulations, elements having similar chemical characteristics (at least, as far as certain properties are concerned), can be compounded into a “super-element. For instance, if the descriptor does not require to distinguish fluorine from chlorine or bromine, these elements can be considered a halogen super-element and accounted for when defining the summation in equation (3). Table 1: The impact of inclusion of atomic charge in the atomic descriptor L on the chemical similarity of local environments. Two template molecules are compared with three test molecules around their protonated nitrogen atoms. The table shows the corresponding local atomic descriptors for the nitrogen atom. The elements not present in the molecules are grouped under the symbol “...” and skipped for clarity. The atomic charge component of the L vector, computed with the OPLS3 force field, is given in the last row of the table.

H ... C N O ... q

Template 1 1.00 0.00 2.53 1.00 0.40 0.00 0.04

Template 2 1.02 0.00 2.50 1.00 0.28 0.00 0.41

Test 1 1.01 0.00 2.52 1.00 0.31 0.00 -0.04

Test 2 1.12 0.00 1.78 1.00 0.28 0.00 0.31

Test 3 1.06 0.00 2.61 1.00 0.27 0.00 -0.22

The second part of descriptor (1), QA , is needed to make the descriptor sensitive to (possibly distant) electronic effects which might influence pKa . For example, a substitution of a hydrogen 10


Page 10 of 50

Page 11 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


atom by a trifluoromethyl group quite far from atom A might lead to a significant change in micro-pKa at A (see, for instance, Ref. 32 ). We also found that certain atomic environments having very different local chemical properties (and micro-pKa ) because of the different formal charge of the central atom might still possess similar geometric surroundings. One such case is the formally uncharged nitrogen atom in a tertiary amine RR0 NR00 and the formally charged nitrogen atom of RR0 N(+) = R000 . Obviously, these nitrogen atoms will have a very different response to protonation (the former is a strong base, and the latter would never be protonated), and yet their local geometric environment can be quite similar (such as three carbon atoms at similar angles and distances in the immediate vicinity). For QA we currently find it sufficient to take an atomic charge at atom A, qA which can be computed with any of the readily available quantum mechanical or even empirical protocols. It should be mentioned that atomic charges have already been used to predict pKa of aminoacids. 33 Since atomic charges are not physical observables it is not possible to define them uniquely. For our purposes the sufficient requirement for the selection of the charge definition is that the method should be able to distinguish between two very different electrostatic environments. See Table 1 and the discussion below for more information on the role of the atomic charge in the local atomic descriptor. A comparison of two local atomic environments A and B usually originating from two different molecules (though in principle A and B are allowed to be located on the same molecule) can be accomplished via computing an overlap LAB between the descriptors LA and LB :

LAB = kLA , LB k h i ~A −M ~ B )2 − wq (qA − qB )2 . = exp −wM (M

(4)

In this formula M is interpreted as a vector, with the difference between two vectors and the square defined in the usual algebraic sense (the latter being a dot product). The “weights” wM and wq are non-negative adjustable parameters which control the contributions of the geometric and electronic components of L to the overlap LAB . The larger the weight the greater the influence

11



Page 12 of 50

of the corresponding component difference on the overlap. By definition, the maximum closeness (overlap) LAB is 1 and the minimum is 0. The similarity of compound X to a whole set of compounds S can be denoted kL(X), L(S)k and is computed as a simple average over all compounds k (NS compounds in total) that constitute set S: kL(X), L(S)k = (1/NS )

NS X

kL(X), L(k)k.

(5)

k=1

In practice the local atomic descriptors for all the compounds in the training set S are precomputed and stored on disk. Table 1 provides an instructive example of the importance of using atomic charges in local atomic descriptors. The charge helps distinguish different atomic environments in situations where the geometric data around a certain atom are very similar, as is the case for all the molecules from Table 1, as far as the protonated nitrogen atom is concerned. Similar geometric environments are characterized by similar geometric (non-charge) components of the LAD vector (see all five columns in Table 1). Consider two “template” molecules from Table 1. These molecules could be part of a training set. We are interested in comparing the environment around the protonated nitrogen atom in three “test” molecules from Table 1 with that of the template molecules. In other words, we are trying to answer the question of which template molecule, 1 or 2, each test molecule is more similar to. When comparing similarity, we are only concerned about the similarity in the vicinity of the protonated nitrogen atom. This appears to be a natural concern when computing local properties driven by local phenomena such as atom protonation/deprotonation. Table 1 shows the similarity measures (overlap) between the local atomic descriptors of the test and the template molecules in two scenarios: when we exclude the charge from the LAD (in such a case assume that wq = 0), and when we do not. Let us first discuss the situation when we ignore the charge while computing the overlap (the numbers in brown). Test 1 is then very close and almost equally close to template 1 and template 2, with a slight preference for template 2. Test 2 is almost equally far from both the templates, but is still a little closer to template 2. And test 3 is very close and almost equally close to both 12


Page 13 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the templates, but is slightly closer to template 2. It is not very reasonable from the chemist’s point of view that test 1 should be almost equally close to both the templates, and above all, that it should lean toward template 2. The protonated nitrogen atom in test 1 belongs to the aromatic system, just as in template 1, but not in template 2. More than that, test 1 and template 1 are very close structurally in other respects, and many chemists might consider them the same group of compounds. Therefore it would be expected that the successful local environment descriptor should place test 1 much closer to template 1 than to template 2. So the LAD without charge fails at this task. When it comes to test 2, the chemist would expect it to be more similar to template 2. In this regard, the chargeless LAD succeeds for test 2, although it does not separate the preference of test 2 for template 2 by much. For test 3, the failure of the chargeless LAD is of the same type as in the situation for test 1. To conclude, the chargeless LAD is ineffective at telling apart molecules with similar geometric environments. Next, we discuss the situation when the charge is taken into account when computing the overlap (the numbers in green). The charges on the protonated nitrogen atoms in templates 1 and 2 are very different (see Table 1). The charges are also very different on the nitrogen atoms in all three tests. This observation makes us hopeful that the charge will be the critical piece of information that will help us more meaningfully classify the test molecules by their preference for the templates. Indeed, taking the charge into consideration resolves almost all the problems observed when dealing with the chargeless LAD. The numbers in green show that test 1 is much more similar to template 1, as it should be, from the chemist’s point of view. Test 2 does not unfortunately show a strong similarity with template 2, but the relative similarity to template 2 is improved. Test 3 is now also much more similar to template 1 than to template 2. The inclusion of charge does not lift all the issues though, as out of tests 1-3 it is test 1 that is still closest to template 2, not test 2, as should be expected. However, it is a big improvement in comparison with the situation when the charge is not in the picture. In our tests of pKa prediction we found that the single atomic charge on the atom of interest is sufficient for LAD to be an effective descriptor of the local environment. We use the charge generation model from the OPLS3 force field 34 in all our calculations. The OPLS3

13



charges are based on CM1A (class IV) charges 35 and use bond correction terms to reproduce the electrostatic potential computed at the Hartree-Fock/6-31G* level of theory as well as solvation free energies on a training set. DFT charges are readily available in our calculations, and could be used too, but we found OPLS3 charges sufficient for the purposes of the present application. The computationally less expensive OPLS3 charges were particularly convenient for developing parameterization in the early versions of our workflow. Should more precise comparison of local environments be needed, additional information about the electrostatics in the vicinity of the atom can be included in the LAD.

3 SHELL MODEL In this section we give a more formal account of the shell model that is used to select parameters in our pKa prediction model. These parameters are tailored to the local environment of the dissociating proton in a chemical compound under study, and will be thus connected to the local atomic descriptor formalism delineated in the preceding section. The previously described Jaguar pKa workflow 1,5,36 predicts pKa values of small molecules using DFT energy calculations in the gas and solution phases. In this workflow, the initially computed so-called “raw” pKa is usually inaccurate and is not directly comparable to experimental pKa . The unavoidable inaccuracies of the DFT gas phase calculations and implicit solvation model which form part of the protocol are present to some degree in all DFT-based pKa prediction methods. They must be modified by empirical corrections dependent on the chemical nature of the local atomic environment of the dissociating proton (functional group). 1 These empirical corrections convert the raw pKa value to the final value that usually exhibits an error of less than 1 pKa unit across diverse chemical structures, especially if conformational effects and the contributions from multiple protonation sites are properly taken into account. Linear regressions or corrections that adjust preliminary pKa values have been used in other pKa prediction protocols. 37–39 The selection of the set of empirical corrections dependent on the dissociating functional group

14


Page 14 of 50

Page 15 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2: A schematic representation of the shell model used in the Jaguar pKa program. Shells are hierarchically organized training sets, each associated with a SMARTS pattern that can be matched by a local environment of a given molecule (shown in green). Shells of lower tiers are considered more general and shells of higher tiers are viewed as more specific. Each shell also has a training set associated with it, as well as the RMSD of the fitting of the training set data to experimental pKa . See text for details. f relies on matching f with the functional group g that is a characteristic of a certain training set S(g). In our protocol we use a hierarchy of overlapping training sets with progressively more general functional groups g, so that often f is matched by several patterns g each associated with its training set S(g). See Figure 2 which illustrates this organization as well as our previous publication for the complete definition of the currently used hierarchical model. 5 The SMARTS patterns of the lower tier shells are progressively more general and “include” the SMARTS patterns of the hierarchically connected higher tier shells. The training sets associated with the lower tier shells include all the training sets of the hierarchically connected higher tier shells, and may have additional training data points not included in these higher tier shells. To return to the previously used example, if f is a pattern described by the chemical formula RN = C(NHR0 )(NHR00 ), and thought of as a guanidine protonatable at the sp2 -hybridized nitrogen atom, it can match the training set of guanidines as well as more general training sets S(g) such as amidines, imines, or even those characterized by the generic nitrogen atom. As a result we will have N matches

15



Page 16 of 50

S(g1 ), S(g2 ), ..., S(gN ), each associated with the set of empirical parameters P (g1 ), P (g2 ), ..., P (gN ). In Figure 2 the shells associated with the matches S(g1 ), S(g2 ), ..., S(gN ) are 1A , 2B , 3E , and 4C , and are marked in green. An application of these sets of parameters leads to N predicted pKa values. In many cases it is not difficult to pick a preferred single pKa value from this set. It can be a value that corresponds to the training set S(gk ) that forms a tighter connection with the group f , the tightness of the connection defined by a suitable automated method, for example, by the number of atoms in the patterns f and gk that overlap. In some cases, however, as demonstrated in Section 1, it is not straightforward to select the best match gk , and therefore to settle on the final single pKa prediction. The raw pKa mentioned earlier in this section is converted to a final pKa value according to model (pattern) gk , as follows:

pKa (gk ) = a(gk )pKa (raw) + b(gk ),

(6)

where a(gk ) and b(gk ) are the only parameters that make up the set P (gk ). For N pattern matches S(g1 ), S(g2 ), ..., S(gN ), as described in the previous paragraph, we will have N corresponding pKa predictions pKa (g1 ), pKa (g2 ), etc. As there is sometimes little chemical or physical justification for preferring one of these predictions over another, a multitude of answers causes increased uncertainty. In practice it is possible to avoid confusion by automatically returning only one value pKa (gk ), picked out of N values. The criterion for this selection can be the number of matched atoms or the root-mean-square deviation (RMSD) of the underlying linear fit of the raw pKa to the experimental data. The idea behind the latter approach, which has been used recently in the Jaguar pKa program, and which we dub the “best shell” method throughout this paper, is this. The shell with the lowest amount of disorder in the corresponding linear fitting should belong to the training set with the greater specificity, so that this shell should be chosen as the most specific among all the shells that are returned as a match. In most of our tests and applications this approach worked very well, 5,9 as the pKa predictions generated by the first shell were the most accurate among all

16


Page 17 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


that matched, when compared to the experimental data. Table 2: Predicted pKa values for a test set of guanidine structures extracted from Ref. 40 The nitrogen atom responsible for the protonation/deprotonation equilibrium is marked in brown. The DFT-based pKa calculations were performed with averaging over 10 conformations with the quantum chemical optimization in the solution phase, according to the method described in Ref. 5 The experimental pKa data refer to measurements in water at room temperature. Outliers are defined as predictions deviating by at least 1.0 pKa unit from the experiment. R1

CH3

R2 NH R4 N

N H

R3

Index 1a 1b 2j 2o 2p 3 4 6 7

R1 Cl Cl H H H Cl H OCH3 Cl

R2 R3 H H H H H H H H H H Cl H H OCH3 Cl H H Cl MUE MSE Outliers

R4 H CH2 CHF2 CH3 CH2 CHF2 CH2 CF3 CH2 CHF2 CH2 CHF2 CH2 CHF2 CH2 CHF2

Best shell 8.72 7.69 10.36 8.75 7.85 7.06 9.08 7.32 6.06 1.20 -1.20 7

Weighted pKa 9.37 8.63 10.57 9.43 8.73 8.17 9.69 8.33 7.25 0.39 -0.39 0

Exp. pKa 9.9 8.9 10.6 9.7 9.2 8.5 10.2 8.9 7.8

However, in some cases picking one value pKa (gk ) from the set {pKa (gk )}k=1,N on the grounds of this criterion (and, possibly, any other similar criterion as long as we pick a single shell) can lead to significant errors and uncertainty in pKa prediction. As an example of such case consider a test set of nine guanidines from our previous work. 5 The molecules from that test set were matched by the shells that described guanidines, amidines, and partly substituted amidines. The shell of guanidines and its associated pKa predictions would be normally the first choice because that shell had the lowest RMSD of data fitting among all these shells and was thus viewed the most specific of the three. The shell of guanidines, however, produced the lowest

17



quality pKa predictions, with mean unsigned errors (MUE) equal to or greater than 1.0 pKa unit, whereas the other two shells, having more dispersion in the training data were far more accurate, giving errors below 0.5 pKa unit. See Table 2 for the details on the guanidine shell predictions. This presumably happened because the shell of guanidines was trained on guanidines of somewhat different chemical composition (mostly aliphatic guanidines) than those from the test set (aromatic guanidines). The molecules from the training sets (shells) of amidines and partly substituted amidines, despite a greater scatter of points in the training data, had more chemical affinity to the test structures than molecules from the guanidines training set (as judged subjectively by the chemist’s eye, and, as we show later in the text, by the local atomic descriptor). Another potential problem with picking a single shell out of N shells is that there can be little preference of one shell over others (because, for example, two or more matched shells can have similar RMSD of data fitting). When such shells have different associated parameters [as used in equation (6)], there is little justification of picking just one shell, and this can mean a substantial uncertainty in pKa prediction. It would be preferable to take a weighed linear combination of the set {pKa (gk )}k=1...N , where weights are evaluated according to the closeness of the test molecule to those from the underlying matched training sets. A scheme proposed in the following section can serve this purpose.

4 WEIGHTED PKA Here we explain how the general idea of local atomic descriptors can be used to generate a weighted pKa prediction. This construction will be free of the problems associated with picking a single “best” shell and its parameters, as described at the end of the previous section. For simplicity of the exposition of the theoretical details we will assume that we are dealing with one protonated and one deprotonated conformation. The methodology that will be constructed for weighted pKa calculation can be trivially extended to multiple conformations (see Section 5). When multiple conformations are used in applications we specify the conformational details explicitly.

18


Page 18 of 50

Page 19 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Best shell method (previous work)

Weighted pKa method (this work)

Figure 3: The conceptual difference between the previously published (best shell) approach and the one presented in the paper. Previously only one matched shell, with the lowest RMSD of the linear fitting to experimental data, would be chosen to select parameters for the final pKa prediction. In the new, improved scheme all matched shells can contribute, in proportion to their chemical closeness to the molecule under investigation X and the quality of their linear fit (value of RMSD). Notice that chemical closeness does not necessarily correlate with RMSD.

19



Page 20 of 50

Suppose we have a set of pKa values {pKa,i }i=1..N associated with the corresponding LADs {Li }i=1..N . Then we can approximate an unknown pKa value associated with an environment of interest characterized by the descriptor Lx as PN pKa,x =

pKa,i kLx , Li k , PN kL , L k x i i

i

(7)

assuming that at least some of the “basis” descriptors Li have a sufficiently large overlap with Lx . Formula (7) represents essentially a weighted interpolation scheme. If both adjustable parameters wM and wq are taken to be zero then kLx , Li k = 1, and pKa,x is a simple average on the set of the pKa,i values (maximum “mixing” of contributions from different values). The larger the values of wM and wq become the less mixing is introduced, with the more and more dominant contributions from those basis descriptors Li that have the greatest overlap with Lx . The optimal weighting parameters wM and wq are chosen after considering a training set of molecules, as described later in the text. It is conceivable to use such a weighted interpolation scheme for other purposes, for example to adjust atomic radii associated with certain local chemical environments in implicit solvation models. Obviously, if all overlaps kLx , Li k in formula (7) are identically zero the expression becomes indeterminate. This situation can occur if we are trying to compute pKa of a molecule in which an atom of an exotic element is being either protonated or deprotonated, and that element is not represented in the training sets. If an overlap of all the basis descriptors with Lx is small then this means that we are attempting to describe a property of an atomic environment through properties of very different atomic environments, a procedure that is going to produce a result of dubious quality. For example, it appears reasonable to be able to approximate the pKa of a secondary amine if we know the pKa ’s of a primary and a tertiary amine. But it would be impolitic to use the latter information to try to predict the pKa of pyridine. Scenarios described in the previous two paragraphs are rare in practice. They can be detected

20


Page 21 of 50

1.3 1.2

Best shell method Weighted pKa method

1.1 Mean Unsigned Error (MUE), pKa units

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0

ls ds ds ds es es es es es es es es es es es es es es es ls no aci aci aci ycl im zin zin zin zol din ilin on lin ilin fon zol id lon oho e c x a a a a i l a m n o ph ino uric ylic ero o ydr er yr pyr pyr y an idi uin y an su etr ioa opo alc p p t th tr t i r r t ox e q m h m i i a a p a rb rb h d yr ti ry r s n a p o ia b ca iou te t c r r se te va

Figure 4: The performance of the best shell (old approach) and weighted pKa (new approach) methods on twenty test sets of small, mainly rigid molecules. The purpose of these tests was to show that the new approach does not deteriorate the results for unproblematic molecules, for which the first shell in the old approach already gives the most accurate pKa prediction.

21



Page 22 of 50

though, at which point it might be judicious to fall back on the parameterization used for the most generic shell. That shell matches every structure and provides “safety in numbers”, in that its parameterization absorbs the behavior of all functional groups and atom types (currently comprising over one thousand data points), and is less likely to lead to a bad prediction than a prediction based on a few very uncertain data points or no data points at all. Using this shell we end up with a single pKa value, in line the idea behind summation (7). Formulas (3)-(5), and (7) can be directly utilized for the purposes of computing weighted pKa once the “basis” of pKa values {pKa (gk )}k=1,N is available. All the shells that are matched against the molecule under investigation [green shells in Figure (2)] are included in the weighting averaging formula (7). Because the shell model is constructed hierarchically, each data point from our training data may belong to more than one of these shells (thus, some shells may overlap with or be subshells of other shells). When we used the above-mentioned formulas directly on our training and test sets with the selection of parameters given below we obtained mostly very accurate pKa predictions, compared to the experiment. However, we noticed that, for a very few functional groups like alcohols, the generic shells which contained a great number of data points and which are known for their little specificity or quality of pKa prediction, had a bit too much closeness to the target structure, according to measure (5). When analyzing the data we realized that it would be beneficial to explicitly include the quality of the training data (that is, RMSD of the linear fit) in our pKa prediction scheme. Indeed, it is not only the closeness of the chemical structures that should define the influence of the training set on the predicted pKa , but also the quality of the training set itself. This idea would allow us to conveniently regulate the amount of contribution from higher- and lower- quality shells to the final value, in particular to mitigate the contribution of overly generic shells with large RMSD’s. The final working formula that we used to compute the pKa of compound X was PN pKa (X) =

k

pKa (gk )kL(X), L(Sk )k e−RMSDk , PN −RMSDk k kL(X), L(Sk )k e

22


(8)

Page 23 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where RMSDk is the RMSD of the linear data fitting of shell Sk . We compute the overlaps kL(X), L(Sk )k of compound X with a set Sk using formula (5). As for X, we always use its protonated form, and always compute the local similarities for the compound X and for the set S at the corresponding protonated heavy atom. It turned out to be unnecessary to introduce an additional parameter in the expression exp (−RMSDk ) in formula (8) because we could achieve very satisfactory results across the whole range of functional groups and test cases by optimizing only four parameters: k1 , k2 , wM , and wq . The new method of pKa calculation based on formula (8) will be referred to as “weighted pKa ” method throughout this paper. A schematic illustration of the principles behind the best shell and weighted pKa methods is given in Figure 3. Aside from the set of scaled covalent atomic radii {RA } for all the elements of the periodic table that we borrowed directly from Grimme and co-authors, 31 there are four parameters in our current model. The parameters k1 and k2 from formula (3) control the influence of surrounding atoms B on atom A. These two parameters have to be optimally adjusted to suppress the influence of faraway atoms, keeping the descriptor local, and at the same time absorb the influence of nearby atoms, giving the descriptors the possibility to capture the identity of the local environment. The parameters wM and wq from formula (4) coordinate the contributions of the geometric factor M defined by formula (2) and the atomic charge q. We need to find balanced values of wM and wq such that we neither overemphasize nor neglect one of the two factors (geometric and electrostatic) in the description of the atomic environment. While seeking the optimal parameters we found that k1 = 16.0 chosen by Grimme and co-authors was already good for our purposes, and we left it at that value. In order to determine the optimal values of the non-linear parameters k2 , wM , and wq we utilized a training set that consisted of 1314 molecules spanning 118 functional group types (shells). Note that the training set consisted of molecules with one relevant functional group only. For these molecules there was no need to combine the computed micro-pKa ’s into a macro-pKa , as explained at the beginning of Section 6. The non-linear optimization processes that we used did not need to be exhaustive as we found

23



Page 24 of 50

that small variations of parameters in the vicinity of the chosen values produced similar results. In the non-linear optimization we first fixed the ratio of wM to wq , while finding the optimal k2 . Even though the individual values of wM and wq influence the overlap (4) only their ratio affects nomalized expressions like (8) which determine the final pKa . In the original work of Grimme and co-workers the k2 value was 4/3. We tested the following values of k2 : 2/3, 1, 4/3, 5/3, 2, and 7/3 in the vicinity of the original value. This produced six overall root mean square errors (RMSE) on the training set: R1 to R6. We then fitted the points (2/3, R1), (1, R2), (4/3, R3), (5/3, R4), (2, R5), and (7/3, R6) to a quadratic function f (x) = ax2 + bx + c, which had a minimum near 5/3. Fixing k2 to be 5/3 we optimized the ratio between wM and wq . In a procedure similar to the one described above, we tested wM and wq to be (2, 5), (3, 5), (7, 10), (4, 5), (9, 10), and (5, 5), which again produced a set of RMSE’s from R1 to R6. An optimization of the associated quadratic function led to the 4/5 as the ratio of wM to wq . Thus, the final optimal parameters are k2 = 5/3, wM = 4.0 and wq = 5.0. All results reported in this work, including those in Table 1, were obtained with these values. We tested our scheme on twenty standard test sets totaling 165 small and medium-sized molecules of moderate flexibility, each of which covers a functional group. These molecules had been used to test Jaguar pKa parameterization in the past. Although the previous Jaguar pKa parameterization schemes had already been performing well on these test sets and had not displayed the sort of problems that led to the development of the novel scheme, it was important to verify that the new approach characterized by equation (8) was not detrimental to the quality of the unproblematic results. Indeed, the test sets showed a very satisfactory performance of the new approach compared to the previous one (see Figure 4 and the Supporting Information). The errors for most functional groups essentially stay constant as we switch to the weighted pKa approach. There are only four functional groups exhibiting a noteworthy change of mean unsigned error of 0.1 pKa unit or more. The weighted pKa method increases the error from 0.33 to 0.51 pKa units for one such group, pyrazines. This increase is acceptable because the final average error is still very modest, consistent with a very accurate pKa prediction. An additional point in favor of the insignificance of this error

24


Page 25 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


increase is that the piperazines test set consisted of only four molecules. The other three functional groups, phenols, tropolones, and alcohols, with 10, 5, and 3 molecules in the test set, respectively, showed a decrease of error by 0.15, 0.46, and 0.20 pKa units, respectively. It is interesting that all the functional groups that showed this slight improvement contained the hydroxy-group.

# --1: 2: 3: 4: 5:

Shell pKa RMSD ------------------------------------------ ----- ----guanidine 6.06 0.320 amidine 7.40 0.662 partly N-substituted amidine 7.88 0.746 Protonation of sp2-like aliphatic nitrogen 6.12 1.751 Protonation of generic atom 6.49 2.264

Similarity ---------14.09% 24.75% 43.72% 10.24% 7.19%

Contribution -----------22.12% 27.60% 44.82% 3.84% 1.62%

Final weighted pKa value is 7.26

# --1: 2: 3: 4:

Shell pKa RMSD ---------------------------------------------- ----- ----phenol excluding o-nitro- and o-nitrosophenols 9.93 0.654 6-member aromatic enol w/ neighboring C(C=O)R 9.05 1.056 Protonation of generic atom 7.19 2.264 Protonation of sp2-like oxygen 6.98 2.568

Similarity ---------37.89% 39.20% 6.40% 16.51%

Contribution -----------55.86% 38.66% 1.89% 3.59%

Final weighted pKa value is 9.43

Figure 5: Jaguar pKa output for two molecules, using 10 conformations per structure with quantum chemical geometry optimization in the solution phase, as described in Ref. 5 The matched shells k, their individual pKa (k) values, RMSDk , the shells’ average similarity to the structures, and their contribution to the final weighted pKa value are shown. The functional groups responsible for the protonation/deprotonation equilibrium are marked in brown. The experimental pKa (k) of the guanidine and phenol molecules are taken from Refs. 40 and, 41 respectively. The Supporting Information provides the SMARTS patterns used for the definition of the shells that appear in this figure. An illustration of the contribution of various shells in typical cases of pKa prediction using Jaguar pKa is given in Figure 5. Two examples with the detailed data related to each shell are shown: one is a guanidine molecule (which is also included in the data in Table 2) and another one is a phenol. In the guanidine molecule example, the shell with the lowest RMSD of the linear fitting (“guanidine”) actually has only a modest similarity around the protonation site to the target 25



molecule. Note that in our old parameterization approach this shell, despite its low structural affinity to the target molecule, would be preferred and would be chosen as the single source of parameters for the final pKa prediction. That choice would be made solely on the basis of the low RMSD of the training set associated with the guanidine shell. Perhaps unsurprisingly, this incoherence results in a very large pKa prediction error of 1.7 pKa units using the old approach. If we utilize the local atomic descriptor at the protonation site and compute the chemical similarity P of shell Sk as kL(X), L(Sk )k/ k kL(X), L(Sk )k in percent we will see that the similarity of the guanidine is only 14.09%, whereas the similarity of the other two shells, “amidine” and “partly N-substituted amidine” is significantly higher, 24.75% and 43.72%, respectively. These two shells have larger pKa values associated with them (7.40 and 7.88, respectively), and, if their greater similarity to the target molecule could be taken into account, they would drive the overall predicted pKa closer to the experimental value, which is 7.8. An accounting for the greater relevance of these two shells could be achieved by the weighted pKa approach. Observe that the similarity of the overly generic shells “protonation of sp2-like aliphatic nitrogen” and “protonation of generic atom”, as measured by the formula indicated above, is non-negligible (10.24% and 7.19%, respectively). The RMSD values associated with these generic shells are much larger than those of the three preceding shells, and stand for a considerable scatter of data points and therefore a lower quality of parameterization. This factor is accounted for by “penalizing” shells by the corresponding exp (−RMSD) multipliers. This suppresses their overall contribution to the P final pKa (computed for shell Sk as kL(X), L(Sk )k e−RMSDk / k kL(X), L(Sk )k e−RMSDk ), and simultaneously increases the contribution of the guanidine shell, as a “reward” for its superior quality of data with low scatter. The final weighted pKa computed by formula (8), balancing the contributions of the shells, is 7.26. This is a much better agreement with the experimental pKa of 7.8 than the value of 6.06 produced by the old, single-shell approach. Having discussed this guanidine molecule in detail, we can skip over the particulars of the other example from Figure 5, the phenol. The similarity and contributions of the two most specific shells in this example are comparable. The only important difference from the previous example is that

26


Page 26 of 50

Page 27 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the impact of the shells on the weighted pKa can be reversed after accounting for their RMSD. The similarity of shell number two to the target molecule in this example was slightly larger than that of shell number one, and would have led to its higher contribution to the weighted pKa if we had not accounted for the shells’ RMSDs. Once the RMSDs are accounted for, the much lower RMSD of the first shell results in the dominant contribution of that shell to the final pKa . The final weighted pKa , computed to be 9.43, is in much better agreement with the experimental value of 9.1 than the pKa of the shell with the lowest RMSD, 9.93, which is what would have been accepted as the final result in our old approach.

5 WORKFLOW DETAILS After having rationalized and exemplified the ideas that constitute the weighted pKa approach, we would like to summarize the workflow’s organization, and, while doing so, list the associated computational and implementational details. The calculation for micro-pKa (characterizing a single protonation/deprotonation process) begins with a specification of a single 3D structure of the molecule supplied by the user. Any conformer can be taken at this point because a conformational search can be applied to this initial structure at a later stage. Our previous tests showed that the pKa results become virtually insensitive to the initial conformation chosen by the user provided one performs a thorough conformational search. 5 Because the user specifies only one structure, the other component of the dissociation process must be algorithmically constructed by either adding or removing a proton. Once both structures that differ by a single proton become available, it is possible to compare them against a standard set of SMARTS patterns from the shell model, and record all the shells that matched. In this way the particular type of the dissociation process is identified. The functional group identification must be performed at this stage because the identity of the functional group might need to trigger some special settings activated later in the workflow (for example, the need to optimize atomic positions

27



of certain zwitterionic species in an implicit solvent model, as opposed to the gas phase). Then, a conformational search integrated in the workflow can be launched. For this step we utilize the program MacroModel 42 with the OPLS2005 force field. 43–45 A more recent and more accurate version of this force field, OPLS3, 34 has not yet been used for generation of conformers in our pKa prediction workflow. The importance of using a thorough conformational search for the quality of pKa prediction results has already been shown. 5,46 As a result of the conformational search we obtain up to n protonated and up to n deprotonated conformations, where n is typically 5 or 10, a number set by the user before the launch of the calculation. The conformational search can return fewer conformations than requested, for example, if the molecule is not flexible enough and possesses only a few distinct conformations. We will denote the final number of conformers used in subsequent calculations by nprot and ndeprot , for the protonated and deprotonated conformations, respectively. The geometry optimization of these conformations is performed with the quantum chemistry program Jaguar 36,47 using the B3LYP/6-31G* level of theory, either in gas or solution phase. For the solution phase calculations we employ the Poisson-Boltzmann finite (PBF) implicit solvation model. 48–50 The gas phase energies and solvation energies required to construct the free energy of the dissociation process are then computed with the B3LYP functional and a mixed cc-pVTZ+/cc-pVTZ basis set. The details of the construction of the free energy ∆G from the energy components and the corresponding thermodynamic cycle are given in earlier publications. 1,5 The energies are computed for all the conformations that are considered, and thus up to n2 ∆Gij values are generated, one for each pair of conformations. Those are converted into pKa (gk ) for each pair of conformations via formula (6). So at this point we have up to n2 gN pKa values, where gN is the total number of matched shells. These pKa values are subsequently corrected by the configurational entropy factors, as explained in Ref., 5 in a procedure that addresses the problem of a continuum of near-isoenergetic conformations. The corrected values are called nano-pKa ’s and denoted pKijnano (gk ), where i stands for the protonated conformation, j for the deprotonated one, and gk for the shell. Now, this multitude of values must be contracted into a single pKa value,

28


Page 28 of 50

Page 29 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: The pKa values predicted by the old and the new pKa parameterization schemes on a set of drug-like molecules from Ref. 51 The nitrogen atom responsible for the protonation/deprotonation equilibrium is marked in brown. The best shell parameterization scheme is already very accurate on this test set. The goal of the present test was to demonstrate that the weighted pKa approach does not substantially change the performance achieved with the best shell model. All the pKa ’s were computed with a geometry optimization in the solution phase using 10 conformations for the protonated and the deprotonated forms, as described in Ref. 5 F O R6 N

NH2

N N H R2

X

R1 R3 R4

Index 1b 66 67 69 109 89 68 88 111 110 14d 70 1a 14a 112 113

R1 CN CN CN CN CN CN CN CN CN CN CN CN Cl Cl Cl Cl

R2 − R5 R3 = F R2 = F R4 = CF3 R5 = CF3 R3 = F R5 = CF3 R5 = CF3 , R3 = F R4 = CF3 , R3 = F R2 = F, R3 = F R2 = F, R3 = F R2 = F, R3 = F R2 = F, R3 = F MUE MSE

R6 CH3 CH3 CH3 CHF2 CH3 CH3 CH2 F CH2 F CH3 CH3 CH3 CH2 F CH3 CH3 CH3 CH3

X O O O O O O O O O O O O O O S S

29

R5

Best shell 9.63 7.97 7.37 7.72 6.68 7.31 7.24 6.17 6.02 5.59 6.36 4.55 9.72 6.45 9.04 6.21 0.25 0.03


Weighted pKa 9.49 7.96 7.42 7.70 6.78 7.34 7.24 6.27 6.15 5.78 6.47 4.74 9.57 6.56 8.84 6.23 0.27 0.06

Exp. pKa 9.7 8.1 7.4 7.3 7.3 7.0 6.7 6.3 5.9 5.8 5.8 5.1 9.8 6.3 9.0 6.1


Page 30 of 50

which was one of the objectives of this work. We first perform a contraction over the indices i, j:

Kamicro (gk )

=

deprot nX

1 Pnprot

j=1

i=1

1/Kijnano (gk )

,

(9)

as is done in our previous publication. 5 We analyze the initial molecule for the existence of functional groups symmetrical to that under investigation. If such symmetry is detected we shift the Kamicro (gk ) values by the corresponding logarithmic factor, which can be derived from formula (11) if we view the protonation of equivalent functional groups as alternative protonation channels. This procedure is well-known and has been applied before, 36,52,53 and we are planning to discuss it in more detail in the context of automated macro-pKa prediction in a future work. Finally, we should be ready to perform a contraction over the shells gk using formula (8). As has already been noted, that formula was derived for a single conformation X in a sense that the protonated local environment of X is compared to that of the compounds from the shells over which the summation is performed. Extension of the previously defined formulas to nprot conformations is trivial and can be done in a number of ways, which in our experience all led to very similar final results. We chose to Boltzmann-average the contributions of multiple conformations to the similarity measure: Pnprot kL(X), L(Sk )k =

i=1

kL(Xi ), L(Sk )k e−Ei /RT , Pnprot −E /RT i i=1 e

(10)

where Ei is the relative energy of the protonated conformation i. Here we work with the protonated forms to be consistent with the convention mentioned in the previous section. So, if structure X is represented by nprot protonated conformations we should use the Boltzmann averaging (10) before using kL(X), L(Sk )k in the final contraction (8). In this case, we compute the overlaps kL(Xi ), L(Sk )k via formula (5). This procedure yields a single micro-pKa value which is reported to the user. Once the micro-pKa is obtained, depending on the molecule, other micro-pKa ’s might need to be computed repeating the process described above, but for a different protonation/deprotonation 30


Page 31 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


channel. Several micro-pKa ’s are then assembled into a final macro-pKa , which is directly comparable to the experimentally measured value, as explained in detail at the beginning of the next section.

6 RESULTS In this section we apply the new weighted pKa method to a few structure types, and compare the quality of its predictions to that of the previously reported method in which only the lowest RMSD shell was used. 5 In our current implementation we still do not account for possibly multiple microstates automatically. Thus, in all the calculations that follow we begin with modeling micro-pKa ’s. The user must define the corresponding protonation/deprotonation channel manually, by specifying either the hydrogen atom that will come off during dissociation, or, alternatively, a heavy atom that will be protonated. A number of molecules from our test sets contain several functional groups or atoms that might be involved in processes that contribute to an experimentally observed macro-pKa . In such cases, we have N prot protonated and N deprot deprotonated microstates, and there are N prot N deprot channels to consider. We set up calculations for each of these protonation/deprotonation channels independently and, having obtained the corresponding micro-pKa ’s, put them together to compute macro-pKa using the formula Kamacro

=

deprot NX

J=1

1 PN prot I=1

micro 1/KIJ

(11)

micro or its equivalents. 2,5,9 In this formula, the dissociation constant KIJ characterizes the dissociation

of protonated microstate I into deprotonated microstate J. In practice, micro-pKa ’s that differ from the dominant micro-pKa by 1.0 pKa unit or more can be ignored, because they would contribute negligibly to the final macro-pKa . In cases of symmetric groups or atoms only a channel for one such symmetric entity needs to be processed explicitly, as the other symmetry-related channels would produce the same micro-pKa ’s. Our implementation actually accounts for such symmetric channels automatically. When no symmetry is present, however, there is currently no alternative for manually setting up each channel and 31



Page 32 of 50

manually assembling the macro-pKa . In our tests N prot and/or N deprot rarely exceed 2, so a manual assembling of macro-pKa is not difficult. However, bookkeeping would become a major problem if the number of microstates were high. An automated generation of possible protonated and deprotonated microstates (or tautomers), along with assembling the results into macro-pKa , is a work in progress and will be presented elsewhere. Among the tests presented below, two types of limiting cases are of particular interest. First, it is desirable and important to prove that the new approach does not deteriorate the already good results obtained with the old approach. Figure 4 has given such a proof for a variety of small, rigid molecules. A very attractive advantage of our previously described method 5 is its ability to accurately predict pKa of large, drug-like molecules owing to a combination of the extensive conformational sampling, accurate DFT energetics, and a sophisticated hierarchical parameterization. 9 This is in contrast to DFT-based pKa prediction methods that are typically applied to small molecules with a single functional group. 54–57 Our conclusion about the method’s ability to accurately predict pKa of more complex molecules is based on tests on dozens of such molecules. Although we observed difficulty with predicting sequential dissociation constants of certain particularly challenging zwitterionic species, there is a strong and promising signal on simpler yet still sufficiently complex molecules. 5,9 It is therefore very important to demonstrate that this advantage is not dissipated as we switch to the even more advanced parameterization scheme proposed in this work. Table 3 shows an application of the old and new parameterization schemes to sixteen drug-like amidines which were used as a test set in our earlier work. 5 The performance of the new weighted pKa method is essentially the same as that of the old method. The mean unsigned error increases only by 0.02 pKa units in comparison with the old approach, which we regard as statistically insignificant, still delivering an impressive 0.27 pKa unit error with zero outliers. Now we move on to the other limiting case. This is where the old method, using only one shell for the final pKa prediction, produces significantly larger errors when the “best” shell’s parameterization is used to compute the final pKa . Other matched shells’ parameterization might

32


Page 33 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


lead to a lower error in the pKa prediction, if their parameters were used, but these shells are disregarded because they have more scatter in their linear fitting (larger RMSDs). This is a situation where the mixed shell approach proposed in this work would be beneficial. Figure 5 and its discussion in the previous section have already presented a demonstration of the benefit of the new approach. Table 2 gives more evidence. The test set presented in the table was problematic for the pKa prediction using our previous approach, and was the only serious issue with the best shell method described in our previous publication. 5 This test set was responsible for unacceptably large individual and mean unsigned pKa prediction errors (over 1.0 pKa unit). The new, weighted pKa , approach fixes this problem. The table demonstrates that the new approach is able to reduce the previously very large mean unsigned error of 1.20 to a mere 0.39. The reason the best shell approach fails for this test set is that the shell with the smallest RMSD of the linear fitting was parameterized largely on aliphatic guanidines. This shell apparently did not have a lot of chemical affinity to the aromatic guanidines from Table 2 (in which the guanidine functional group was linked to an aromatic ring). Impressively, the LAD sensed a low chemical similarity of the structures from the lowest RMSD shell. Consequently, the weighed pKa method readjusted contributions of the shells, assigning only little weight to the first shell, and giving most weight to the shell of “partly N-substituted guanidines”. The latter, most relevant shell, comes only third in RMSD ranking. See Figure 5 that gives the details on shell contributions for one structure from Table 2 (which is the structure with index 7) . That shell’s parameters would not be used for the final pKa calculation in our previous approach, and it would be difficult to prefer without either a preliminary knowledge of the experimental results or a careful inspection of the training sets by eye.

7

COMPARISON WITH OTHER METHODS

The main purpose of our work was a development of a superior parameterization technique based on a novel local descriptor. We have shown how this new scheme conceptually simplifies the output

33



and interpretation of the generated pKa data and in certain cases greatly improves the accuracy of results. Competition for higher accuracy with other types of pKa prediction methods was not a primary objective of our work. Nevertheless, it would be useful to put the currently proposed method in the context of other approaches by comparing the accuracy of different methods on test sets. Such a step should provide an additional validation of our approach, on top of what was already described in the previous sections. It is not always easy to decide what constitutes a fair comparison, as different methods are often designed for different goals and might have different limitations. A given test set might expose these limitations seemingly to a disadvantage of some methods. When interpreting comparisons of methods’ performance it is also useful to keep in mind that the new test sets and results accompanying novel methods presented in the literature (doubtless including those from our own works) unavoidably suffer from so-called publication bias. 58 That is to say, the test sets and (the best) results accompanying the novel method tend to be chosen, almost by necessity, so as to present the method in a positive light. This is significant because such test sets, if used for testing a competing approach, are essentially optimized to give advantage to the associated published methods that one competes against. This problem occurs when methods of interest are not easily available for testing on a new, presumably less biased test set. We begin our comparisons with other pKa methods with Figure 6. We declare that the test set presented in it, in view of what was said in the preceding paragraph, is of our own composition. In this application, we studied four relatively rigid but still sufficiently large morphine-like structures. Their experimental pKa ’s are very sensitive to the stereochemistry of the hydroxyl-group which neighbors the amino-group responsible for the protonation/deprotonation equilibrium. We applied four different computational methods to predict the pKa of these molecules: both the best shell and the weighted pKa methods used in Jaguar pKa, as well as empirical approaches incorporated into the programs Epik 7,59,60 developed by Schrödinger and the pKa module of the interface Chemicalize 61 developed by ChemAxon. All four methods delivered impressive pKa predictions with low mean unsigned errors of 0.55, 0.58, 0.54, and 0.72 pKa units, respectively. The very small

34


Page 34 of 50

Page 35 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


differences in the first three errors should be viewed as statistically insignificant. The weighted pKa approach gives virtually the same quality result as the best shell approach. Chemicalize is slightly less accurate on this test set than Epik, but in view of the equally incorrect trend that both empirical approaches show in Figure 6, the larger error of Chemicalize appears fortuitous. Importantly, even though the average errors of the DFT-based and empirical methods are quite similar and small, the trends in pKa change with the alteration of the amino group and the varying stereochemistry of the neighboring hydroxyl are impressively reproduced only by the DFT-based approach. Both the empirical methods are actually completely insensitive to the stereochemistry of the hydroxy-group. Empirical methods are expected to give low average pKa prediction errors on uncomplicated, relatively rigid amines, because the amino-group has been extensively studied and parameterized. However, it is the trends in pKa variation that are also important in examples like the test set from Figure 6, and only the DFT-based approach describes both the trend and the individual pKa values accurately. Our second test for the purpose of comparison with other methods comes from the recent publication by Haworth and co-workers. 46 It is a set of primary, secondary, and tertiary amines comprising 44 molecules and 53 individual pKa ’s (as consecutive dissociation constants are considered for diamines). We compare the results of our method on this test set with the results of what appears to be one of the most accurate methods of Haworth and co-workers called PEX HL-TC, which employs the highly accurate G3 theory 63 for free energy calculations and a sophisticated coverage of the conformational space. Even though certain molecules from the test set from Ref. 46 are large and flexible they are not as large and flexible as some other molecules studied in our work (see in particular our Table 3). A number of molecules from this set are small and a few of them turned out to have been included in our training sets for developing empirical parameters used in our approach. Since the overlap with our training set is small and since the empirical parameters in our approach are only some of the several factors influencing the final predicted pKa , these molecules are not going to substantially sway the overall result of the test. To remove all doubt, we have also performed a comparison between our method and PEX HL-TC on a subset

35



Figure 6: Application of computational pKa prediction models to a set of heterocyclic structures in which stereochemistry has a strong impact on experimental pKa . The nitrogen atoms responsible for the protonation/deprotonation equilibrium are marked in brown. The DFT-based pKa calculations were performed with averaging over 10 conformations with the quantum chemical optimization in the solution phase, according to the method described in Ref. 5 The experimental pKa data were taken from Ref. 62 See text for discussion.

36


Page 36 of 50

Page 37 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of 34 structures that omitted the molecules used in our training. The MUE of the PEX HL-TC for the full set comprising 44 species is 0.38 pKa units. Our weighted pKa approach with 10 conformations and geometry optimization in solution phase, applied to the full test set, gives 0.47 pKa units for the mean unsigned error, which can be considered a good performance. The reduced test set of 34 structures, defined above, gives an almost identical result: 0.37 and 0.47 units for the mean unsigned error. The detailed results on each structure from the test set of Haworth and co-workers are provided in the Supporting Information. Table 4: Comparison of performance of various pKa prediction methods on a test set of Yu and co-workers 64 which is composed of 1143 small and medium-sized molecules. See the original work and reference therein for the details on the QC, r-QC, SPARC, and ACD methods. The individual results for each molecule are given in the Supporting Information.

MUE MSE

Best shell 0.50 -0.11

Weighted pKa 0.49 -0.15

QC 0.59 -0.07

r-QC 0.54 0.00

SPARC 0.39 0.00

ACD 0.30 -0.01

The third and last test comparing the performance of our method with other pKa prediction methods is perhaps the most challenging of all those discussed so far, and it comes from a fairly recent work by Yu and co-workers. 64 The test set contains 1143 small and medium-sized molecules (580 oxygen acids and 563 nitrogen bases), and as many pKa data points. Yu and co-workers applied several pKa prediction protocols to this test set and compared their performance. We applied our best shell and weighted pKa methods to this test set and compared the results with those of other methods. The detailed information for every molecule from the test set is given in Supporting Information, while Table 4 summarizes the mean signed and unsigned errors. Our calculations used up to 10 conformations per protonation form with subsequent geometry optimization in gas phase. We took the experimental values from the paper of Yu and co-workers even though slightly different experimental pKa values were found by us in other sources for some data points. The mean unsigned errors produced by our two methods are similar and fall within the range of errors produced by the other methods listed in Table 4 . In particular, for this test set, our methods are slightly more accurate than the semiempirical methods QC and r-QC and are somewhat less 37



accurate than the empirical methods SPARC and ACD (see Ref. 64 and references therein for a definition and discussion of these methods). The mean signed errors indicate that our methods tend do underpredict pKa by values about 0.1 unit, the QC method shows a bit smaller systematic underprediction, whereas the empirical methods show no systematic shift. The slightly negative MSE produced by our methods can be perhaps explained by the construction of training sets of for our shells. The presence of molecules with electron-withdrawing groups in these sets is likely slightly larger than encountered in practical test sets. This is further corroborated by the results on amidines shown in Tables 2 and 3. The first test set containing molecules with moderate electron withdrawing groups shows a strong tendency of our methods to underpredict pKa . In contrast, the amidine test set from Table 3 composed mostly of many molecules with strong electron withdrawing groups and consequently systematically lower experimental pKa ’s shows a slight tendency of our methods toward overprediction. In this section we have shown applications of our novel method to three test sets of varying complexity and compared the results to those of several methods previously published by other groups. The details of the accuracy of the pKa predictions differ from application to application, but the overall conclusion is that the performance of the weighted pKa method presented in this work is comparable to the standards of other modern pKa prediction methods.

8 DISCUSSION AND OUTLOOK The innovations proposed in this paper are twofold. First, in Section 2 we developed a descriptor of local chemical environment that enabled us to evaluate chemical similarity between two local environments (in either the same or different structures). We have shown that it was also possible to compare the local environment of a given structure with that of a set of structures, such as a training set. Second, we improved the previously reported application of empirical parameters to the raw DFT energetics in a process that produces final pKa values comparable to the experiment. The new scheme, described in Section 4, combines parameters from a number of training sets in

38


Page 38 of 50

Page 39 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


proportion to the chemical similarity between the given structure and the structures constituting the training set. Even though that approach produced generally good results, we found that the results could be further improved, and some outliers in the final data can be eliminated, by weighing in the quality of the training sets. As a consequence, the contribution of the training set to the parameters used to compute the final pKa value is influenced by both the structural similarity of the molecules in the training set and the quality of the data fitting on that set. This new scheme is expected to be smoother in a sense that it is less likely to lead to unjustified, abrupt changes in pKa predictions upon modification of the chemical structure, caused by a sudden shift in the choice of empirical parameters. It is more accurate because it eliminates large errors originating from an association of a given structure with a single training set which represents it inadequately (see Table 2). The accurate results produced by the new model are exemplified on a number of tests. By construction, the new model would not be expected to substantially alter the already good results that the previously published approach was capable of delivering, and instead would be expected to fix the situations in which the old model produced large errors. The tests convincingly demonstrate that this is the case. Nevertheless, certain deficiencies in our presently used workflow and the parameterization scheme remain. The first concern is the quality and chemical scope of the training sets in the shell model. If all shells in the shell model have little chemical similarity to the target molecule, then the new scheme, even though it will formally work, is not expected to produce accurate results. A way to address this concern has already been discussed in the text. It presumes falling back on the generic parameterization, which is not an entirely satisfactory solution. In the following paragraph, when discussing one more deficiency of the current workflow, we consider another way of addressing this problem. The construction of the shell model, a hierarchical organization of SMARTS patterns, is at present done by a human researcher. This arrangement is inefficient and largely subjective. There are no clear rules of which SMARTS patterns should define each shell, and how to organize the shells into a hierarchy. The shell model, while making Jaguar pKa applicable to essentially all

39



Page 40 of 50

chemical space with an automatic adjustment of specificity of parameterization depending on the functional group, is difficult to construct, maintain, and expand. The operation of adding new data points to the model is laborious and not entirely algorithmic. While the guidelines for adding new data points exist, it is still up to the human researcher’s final judgment to classify the data point into one or more of the shells. A question arises of whether it is possible to dispense with SMARTS patterns in the construction of the pKa parameterization altogether because one could imagine that matching local environments of molecules could instead be performed with LADs. A training set of molecules with known experimental pKa and precomputed LADs would be stored on disk. The LAD at the molecule of interest would be compared with the LADs of the training set. A number of molecules of sufficient similarity would be selected as a result of this matching, and linear parameters such as a and b from equation (6) would be constructed “on the fly”. Alternatively, all the data points in the large training set could contribute to the parameterization, but in proportion to their chemical similarity. In this arrangement, the parameterization would be tailored to the molecule of interest. This scheme would obviate the need to construct and maintain the hierarchical organization of shells. This is an interesting possibility that deviates radically from the parameterization scheme which already proved to be successful. Because it would require additional careful testing and validation we leave its exploration for future work. The idea to replace pattern matching based on SMARTS patterns with that based on LADs might also be applicable in other areas of computational research. We emphasize that this replacement is expected to work well only in cases of local structure matching. There are models and implementations in which atomic or other local parameters are assigned on the basis of matching against a library of SMARTS patterns. 65–67 If we assume that the use of LADs is efficacious and can replace the SMARTS pattern matching, we would have to deal with a library of 3D molecular structures and not SMARTS patterns. A library of molecules is easier to maintain and eliminates multiple problems associated with SMARTS patterns to which we referred earlier in the text, such as problems with multiple possible resonance structures and problems with the definition of covalent

40


Page 41 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


bonding in nontrivial molecules. This idea still needs to be tested in practice, but it appears that parameterization models based on LAD hold more promise than a single application to the calculation of pKa presented in this work.

9 ACKNOWLEDGMENT The authors would like to thank ChemAxon, Ltd. for the permission to publish Chemicalize pKa results.

10 SUPPORTING INFORMATION Detailed pKa data for test sets shown in Figure 4 are provided. The detailed results on the test set of Haworth and co-workers 46 as well as a breakdown of results from Table 4 by molecule are also available.

11 CORRESPONDING AUTHOR Art D. Bochevarov address: Schrödinger, Inc., 120 West 45th St, New York, New York 10036 e-mail: [email protected]

41



References (1) Klicić, J. J.; Friesner, R. A.; Liu, S.-Y.; Guida, W. C. Accurate Prediction of Acidity Constants in Aqueous Solution via Density Functional Theory and Self-Consistent Reaction Field Methods. J. Phys. Chem. A 2002, 106, 1327–1335. (2) Jang, Y. H.; Goddard III, W. A.; Noyes, K. T.; Sowers, L. C.; Hwang, S.; Chung, D. S. pKa Values of Guanine in Water: Density Functional Theory Calculations Combined with PoissonBoltzmann ContinuumSolvation Model. J. Phys. Chem. B 2003, 107, 344–357. (3) Eckert, F.; Klamt, A. Accurate Prediction of Basicity in Aqueous Solution with COSMO-RS. J. Comp. Chem. 2006, 27, 11–19. (4) Riojas, A. G.; Wilson, A. K. Solv-ccCA: Implicit Solvation and the Correlation Consistent Composite Approach for the Determination of pKa. J. Chem. Theory Comput. 2014, 10, 1500–1510. (5) Bochevarov, A. D.; Watson, M. A.; Greenwood, J. R.; Philipp, D. M. Multiconformation, Density Functional Theory-Based pKa Prediction in Application to Large, Flexible Organic Molecules with Diverse Functional Groups. J. Chem. Theory Comput. 2016, 12, 6001–6019. (6) Hansen, N. T.; Kouskoumvekaki, I.; Jørgensen, F. S.; Brunak, S.; Jónsdóttir, S. O. Prediction of pH-Dependent Aqueous Solubility of Druglike Molecules. J. Chem. Inf. Model. 2006, 46, 2601–2609. (7) Shelley, J. C.; Cholleti, A.; Frye, L.; Greenwood, J. R.; Timlin, M. R.; Uchimaya, M. Epik: a Software Program for pK( a ) Prediction and Protonation State Generation for Drug-like Molecules. J. Comput. Aided Mol. Des. 2010, 21, 681–691. (8) Fraczkiewicz, R.; Lobell, M.; Göller, A. H.; Krenz, U.; Schoenneis, R.; Clark, R. D.; Hillisch, A. Best of Both Worlds: Combining Pharma Data and State of the Art Modeling Technology To Improve in Silico pKa Prediction. J. Chem. Inf. Model. 2015, 55, 389–397. 42


Page 42 of 50

Page 43 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(9) Philipp, D. M.; Watson, M. A.; Yu, H. S.; Steinbrecher, T. B.; Bochevarov, A. D. Quantum Chemical pKa Prediction for Complex Organic Molecules. Int. J. Quantum Chem. 2017, e25561. (10) Willett, P.;

Barnard, J. M.;

Downs, G. M. Chemical Similarity Searching.

J. Chem. Inf. Comput. Sci. 1998, 38, 983–996. (11) Streitwieser, A.; Heathcock, C. H. Introduction to Organic Chemistry; MacMillan: New York, 1985. (12) Wade, L. G. Phenol; Encyclopaedia Britannica, 2008. (13) McMurry, J. Fundamentals of Organic Chemistry; Cengage Learning: Belmont, CA, 2010. (14) Lemke, T. L. Review of Organic Functional Groups: Introduction to Medicinal Organic Chemistry; Lippincott Williams & Wilkins: Baltimore, MD, 2003. (15) Vollhardt, P.; Schore, N. Organic Chemistry: Structure and Function; Cengage Learning: New York, 2014. (16) De Proft, F.; Geerlings, P. Conceptual and Computational DFT in the Study of Aromaticity. Chem. Rev. 2001, 101, 1451–1464. (17) Friedrichsen, W.; Traulsen, T.; Elguero, J.; Katritzky, A. R. Advances in Heterocyclic Chemistry; 2000; Vol. 76; pp 85–156. (18) Bartók, A. P.; Kondor, R.; Csányi, G. On Representing Chemical Environments. Phys. Rev. B 2013, 87, 184115. (19) Ferré, G.; Maillet, J.-B.; Stoltz, G. Permutation-invariant Distance Between Atomic Configurations. J. Chem. Phys. 2015, 143, 104114.

43



(20) von Lilienfeld, O. A.; Ramakrishnan, R.; Rupp, M.; Knoll, A. Fourier Series of Atomic Radial Distribution Functions: a Molecular Fingerprint for Machine Learning Models of Quantum Chemical Properties. Int. J. Quant. Chem. 2015, 115, 1084–1093. (21) Khorshidi, A.; Peterson, A. A. Amp: A Modular Approach to Machine Learning in Atomistic Simulations. Comp. Phys. Comm. 2016, 207, 310–324. (22) Nohair, M.; Zakarya, D.; Berrada, A. Autocorrelation Method Adapted To Generate New Atomic Environments: Application for the Prediction of 13-C Chemical Shifts of Alkanes. J. Chem. Inf. Comput. Sci. 2002, 42, 586–591. (23) Smith, J. S.; Isayev, O.; Roitberg, A. E. ANI-1: an Extensible Neural Network Potential with DFT Accuracy at Force Field Computational Cost. Chem. Sci. 2017, 8, 3192–3203. (24) James, C. A.; Weininger, D. Daylight Theory Manual, Version 4.9. Daylight Chemical Information Systems, Inc.: Laguna Niguel, CA, 2011. (25) Bietz, S.; Schomburg, K. T.; Hilbig, M.; Rarey, M. Discriminative Chemical Patterns: Automatic and Interactive Design. J. Chem. Inf. Model. 2015, 55, 1535–1546. (26) Hähnke, V. D.; Bolton, E. E.; Bryant, S. H. PubChem Atom Environments. J. Cheminform. 2015, 7, 41. (27) Lee, A. C.; Yu, J.-Y.; Crippen, M. pKa Prediction of Monoprotic Small Molecules the SMARTS Way. J. Chem. Inf. Model. 2008, 48, 2042–2053. (28) Guba, W.; Meyder, A.; Rarey, M.; Hert, J. Torsion Library Reloaded: A New Version of Expert-Derived SMARTS Rules for Assessing Conformations of Small Molecules. J. Chem. Inf. Model. 2016, 56, 1–5. (29) Proschak, E.; Wegner, J. K.; Schüller, A.; Schneider, G.; Feschner, U. Molecular Query Language (MQL)A Context-Free Grammar for Substructure Matching. J. Chem. Inf. Model. 2007, 47, 295–301. 44


Page 44 of 50

Page 45 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(30) Homer, R. W.; Swanson, J.; Jilek, R. J.; Hurst, T.; Clark, R. D. SYBYL Line Notation (SLN): A Single Notation To Represent Chemical Structures, Queries, Reactions, and Virtual Libraries. J. Chem. Inf. Model. 2008, 48, 2294–2307. (31) Grimme, S.; Anthony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (32) Baranov, M. S.; Solntsev, K. M.; Lukyanova, K. A.; Yampolsky, I. V. A Synthetic Approach to GFP Chromophore Analogs from 3-Azidocinnamates. Role of Methyl Rotors in Chromophore Photophysics. Chem. Commun. 2013, 49, 5778–5780. (33) Ugur, I.; Marion, A.; Parant, S.; Jensen, J. H.; Monard, G. Rationalization of the pKa Values of Alcohols and Thiols Using Atomic Charge Descriptors and Its Application to the Prediction of Amino Acid pKas. J. Chem. Inf. Model. 2014, 54, 2200–2213. (34) Harder, E.; Damm, W.; Maple, J.; Wu, C.; Reboul, M.; Xiang, J. Y.; Wang, L.; Lupyan, D.; Dahlgren, M. K.; Knight, J. L.; Kaus, J. W.; Cerutti, D. S.; Krilov, G.; Jorgensen, W. L.; Abel, R.; Friesner, R. A. OPLS3: A Force Field Providing Broad Coverage of Drug-like Small Molecules and Proteins. J. Chem. Theor. Comput. 2016, 12, 281–296. (35) Storer, J. W.; Giesen, D. W.; Cramer, C. J.; Truhlar, D. G. Class IV Charge Models: A New Semiempirical Approach in Quantum Chemistry. J. Comp.-Aid. Mol. Des. 1995, 9, 87–110. (36) Bochevarov, A. D.; Harder, E.; Hughes, T. F.; Greenwood, J. R.; Braden, D. A.; Philipp, D. M.; Rinaldo, D.; Halls, M. D.; Zhang, Z.; Friesner, R. A. Jaguar: A High-performance Quantum Chemistry Software Program with Strengths in Life and Materials Sciences. Int. J. Quant. Chem. 2013, 113, 2110–2142. (37) Klamt, A.; Eckert, F.; Diedenhofen, M.; Beck, M. E. First Principles Calculations of Aqueous pKa Values for Organic and Inorganic Acids Using COSMORS Reveal an Inconsistency in the Slope of the pKa Scale. J. Phys. Chem. A 2003, 107, 9380–9386. 45



(38) Brown, T. N.; Mora-Diez, N. Computational Determination of Aqueous pKa Values of Protonated Benzimidazoles (Part 1). J. Phys. Chem. B 2006, 110, 9270–9279. (39) Farràs, P.; Teixidor, F.; Branchadell, V. Prediction of pKa Values of nido-Carboranes by Density Functional Theory Methods. Inorg. Chem. 2006, 45, 7947–7954. (40) Peters, J.-U.; Lübbers, T.; Alanine, A.; Kolczewski, S.; Blasco, F.; Steward, L. Cyclic Guanidines as Dual 5-HT5A/5-HT7 Receptor Ligands: Optimising Brain Penetration. Bioorg. Med. Chem. Lett. 2008, 18, 262–266. (41) Maugh II, T.; Bruice, T. C. Role of Intramolecular Bifunctional Catalysis of Ester Hydrolysis in Water. J. Am. Chem. Soc. 1971, 93, 3237–3248. (42) Schrödinger Release 2017-2: MacroModel. Schrödinger, LLC: New York, NY, 2017. (43) Jorgensen, W. L.; Tirado-Rives, J. The OPLS [Optimized Potentials for Liquid Simulations] Potential Functions for Proteins, Energy Minimizations for Crystals of Cyclic Peptides and Crambin. J. Am. Chem. Soc. 1988, 110, 1657–1666. (44) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225–11236. (45) Shivakumar, D.; Williams, J.; Wu, Y.; Damm, W.; Shelley, J.; Sherman, W. Prediction of Absolute Solvation Free Energies using Molecular Dynamics Free Energy Perturbation and the OPLS Force Field. J. Chem. Theory Comput. 2010, 6, 1509–1519. (46) Haworth, N. L.; Wang, Q.; Coote, M. L. Modeling Flexible Molecules in Solution: A pKa Case Study. J. Phys. Chem. A 2017, 121, 5217–5225. (47) Schrödinger Release 2017-2: Jaguar. Schrödinger, LLC: New York, NY, 2017. (48) Cortis, C. M.; Friesner, R. A. An Automatic 3D Finite Element Mesh Generation System for the Poisson-Boltzman Equation. J. Comp. Chem. 1997, 18, 1570–1590. 46


Page 46 of 50

Page 47 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(49) Cortis, C. M.; Friesner, R. A. Numerical Solution of the Poisson-Boltzmann Equation Using Tetrahedral Finite-Element Meshes. J. Comp. Chem. 1997, 18, 1591–1608. (50) Marten, B.; Kim, K.; Cortis, C.; Friesner, R. A.; Murphy, R. B.; Ringnalda, M. N.; Sitkoff, D.;

Honig, B. New Model for Calculation of Solvation Free Energies:

Correction of Self-Consistent Reaction Field Continuum Dielectric Theory for Short-Range Hydrogen-Bonding Effects. J. Phys. Chem. 1996, 100, 11775–11788. (51) Hilpert, H.;

Guba, W.;

Woltering, T. J.;

Wostl, W.;

Pinard, E.;

Mauser, H.;

Mayweg, A. V.; Rogers-Evans, M.; Humm, R.; Krummenacher, D.; Muser, T.; Schnider, C.; Jacobsen, H.; Ozmen, L.; Bergadano, A.; Banner, D. W.; Hochstrasser, R.; Kuglstatter, A.; David-Pierson, P.; Polara, H. F. A.; Narquizian, R. β-Secretase (BACE1) Inhibitors with High in Vivo Efficacy Suitable for Clinical Evaluation in Alzheimers Disease. J. Med. Chem. 2013, 56, 3980–3995. (52) Khalili, F.; Henni, A.; East, A. L. L. Entropy contributions in pKa computation: Application to alkanolamines and piperazines. J. Mol. Structure: THEOCHEM 2009, 916, 1–9. (53) Jang, Y. H.; Hwang, S.; Chang, S. B.; Ku, J.; Chung, D. S. Acid Dissociation Constants of Melamine Derivatives from Density Functional Theory Calculations. J. Phys. Chem. A 2009, 113, 13036–13040. (54) Casasnovas, R.; Fernández, D.; Ortega-Castro, J.; Frau, J.; Donoso, J.; Muñoz, F. Avoiding Gas-phase Calculations in Theoretical pKa Predictions. Theor. Chem. Acc. 2011, 130, 1–13. (55) Rebollar-Zepeda, A. M.; Galano, A. First Principles Calculations of pKa Values of Amines in Aqueous Solution: Application to Neurotransmitters. Int. J. Quant. Chem. 2012, 112, 3449–3460. (56) Thapa, B.; Schlegel, H. B. Theoretical Calculation of pKa’s of Selenols in Aqueous Solution Using an Implicit Solvation Model and Explicit Water Molecules. J. Phys. Chem. A 2016, 120, 8916–8922. 47



(57) Thapa, B.; Schlegel, H. B. Density Functional Theory Calculation of pKa’s of Thiols in Aqueous Solution Using Explicit Water Molecules and the Polarizable Continuum Model. J. Phys. Chem. A 2016, 120, 5726–5735. (58) Fanelli, D. “Positive” Results Increase Down the Hierarchy of the Sciences. PLoS ONE e10068. (59) Schrödinger Release 2017-2: Epik. Schrödinger, LLC: New York, NY, 2017. (60) Greenwood, J. R.; Calkins, D.; Sullivan, A. P.; Shelley, J. C. Towards the Comprehensive, Rapid, and Accurate Prediction of the Favorable Tautomeric States of Drug-like Molecules in Aqueous Solution. (61) Chemicalize. ChemAxon (http://www.chemaxon.com), 2017. (62) Rapoport, H.; Masamune, S. The Stereochemistry of 10-Hydroxycodeine Derivatives. J. Am. Chem. Soc. 1955, 77, 4330–4334. (63) Henry, D. J.; Sullivan, M. B.; Radom, L. G3-RAD and G3X-RAD: Modified Gaussian-3 (G3) and Gaussian-3X (G3X) Procedures for Radical Thermochemistry. J. Chem. Phys. 2003, 118, 4849–4860. (64) Yu, H.; Kühne, R.; Ebert, R.-U.; Schüu¨ rmann, G. Comparative Analysis of QSAR Models for Predicting pKa of Organic Oxygen Acids and Nitrogen Bases from Molecular Structure. J. Chem. Inf. Model. 2010, 50, 1949–1960. (65) Hou, T. J.; Xu, X. J. ADME Evaluation in Drug Discovery. 2. Prediction of Partition Coefficient by Atom-Additive Approach Based on Atom-Weighted Solvent Accessible Surface Areas. J. Chem. Inf. Comput. Sci. 2003, 43, 1058–1067. (66) Zhang, W.; Hou, T.; Xu, X. New Born Radii Deriving Method for Generalized Born Model. J. Chem. Inf. Model. 2005, 45, 88–93.

48


Page 48 of 50

Page 49 of 50 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(67) Kenny, P. W.; Montanari, C. A.; Prokopczyk, I. M. ClogPalk: a Method for Predicting Alkane/water Partition Coefficient. J. Comp.-Aid. Mol. Des. 2013, 27, 389–402.

49



42x33mm (300 x 300 DPI)


Page 50 of 50

A Weighted Averaging Scheme and a Local Atomic Descriptor for pKa

Recommend Documents