Multiconformation, Density Functional Theory ... - ACS Publications

Subscriber access provided by NEW YORK UNIV

Article

Multiconformation, density functional theorybased pKa prediction in application to large, flexible organic molecules with diverse functional groups Art D. Bochevarov, Mark A. Watson, Jeremy R. Greenwood, and Dean M. Philipp J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00805 • Publication Date (Web): 10 Nov 2016 Downloaded from http://pubs.acs.org on November 29, 2016

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 Theory and Computation 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 24

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 Theory and Computation

Multiconformation, density functional theory-based pKa prediction in application to large, flexible organic molecules with diverse functional groups Art D. Bochevarov,∗ Mark A. Watson, and Jeremy R. Greenwood Schrödinger, Inc., 120 West 45th St, New York, New York 10036

Dean M. Philipp Schrödinger, Inc., 101 SW Main Street, Suite 1300, Portland, Oregon 97204 We consider the conformational flexibility of molecules and its implications for micro- and macro-pKa . The corresponding formulas are derived and discussed against the background of a comprehensive scientific and algorithmic description of the latest version of our computer program Jaguar pKa, a density functional theory-based pKa predictor, which is now capable of acting on multiple conformations explicitly. Jaguar pKa is essentially a complex computational workflow incorporating research and technologies from the fields of cheminformatics, molecular mechanics, quantum mechanics, and implicit solvation models. The workflow also makes use of automatically applied empirical corrections which account for the systematic errors resulting from the neglect of explicit solvent interactions in the algorithm’s implicit solvent model. Applications of our program to large, flexible organic molecules representing several classes of functional groups are shown, with a particular emphasis in illustrations laid on drug-like molecules. It is demonstrated that a combination of aggressive conformational search and an explicit consideration of multiple conformations nearly eliminates the dependence of results on the initially chosen conformation. In certain cases this leads to unprecedented accuracy, which is sufficient for distinguishing stereoisomers that have slightly different pKa values. An application of Jaguar pKa to proton sponges, the pKa of which are strongly influenced by steric effects, showcases the advantages that pKa predictors based on quantum mechanical calculations have over similar empirical programs. Keywords:

I. INTRODUCTION

Despite a thorough general understanding of the physical processes that determine the protonation equilibria of organic molecules in aqueous solutions, a widely-applicable and robust computational prediction of the negative logarithm of the protonation equilibrium constant, pKa , remains a challenge.1–7 The pKa of a chemical compound is an important characteristic in drug discovery8–10 and materials science11–13 applications. Numerous studies have shown that it is possible to accurately compute pKa using approaches based on either quantum mechanics14–18 or cheminformatics.19,20 However, a large number of such studies focus on a narrow set of functional groups or families of chemical structures,14–16,18,20–29 and it is not clear if the accuracy of pKa prediction claimed in these works would persevere when using the same approaches for other types of structures that inevitably come up in practical applications. There is also a concern about using somewhat different computational protocols for predicting pKa of different classes of compounds,14,15 or using different experimental proton solvation free energy values in the models that are meant to generate absolute pKa ,2,14,23,30,31 without the need for empirical adjustments. These issues raise doubts about the availability of a truly universal computational method for pKa prediction. Computational methods that predict pKa can be divided into two broad classes: methods that utilize atomic coordinates (three-dimensional structures) and methods that utilize atom connectivity information (two-dimensional structures). The first class of methods, which we will be calling 3Dmethods, deals with conformations explicitly. The use of 3D structures, and consequently conformations, can be an ad-

vantage and a disadvantage for pKa predictions. An advantage is that such an approach, by virtue of having access to atomic coordinates, can be made responsive to steric, conformational, and other subtle effects such as intramolecular hydrogen bonds or π-π interactions, all of which might have a significant impact on pKa . A disadvantage is that a thorough and accurate accounting for conformations is difficult and computationally expensive to perform, so that a failure to do this part of the workflow properly or a temptation to use a less computationally expensive conformational search (CS) might lead to large pKa errors for conformationallysensitive structures. The second class of methods that deals only with atom connectivities and ignores the conformational issues entirely can be paradoxically more stable when handling conformationally-sensitive structures. This can happen either because in such methods the 3D effects are factored in the underlying extensive parameterization, or because the error, even if large, generated by these methods for the same type of structures is not affected by the essentially random initial conformation, and tends to be more systematic. Even if pKa can be evaluated accurately on sets of conformationally flexible molecules with diverse functional groups, this might be still insufficient for achieving a good agreement with the experimental pKa in more challending cases. Such challenging systems that require additional treatment from the point of view of the pKa predictor are: molecules containing multiple functional groups which can be protonated or deprotonated, molecules that exist in different tautomeric forms (that is, forms that differ by the localization of a single hydrogen atom or proton), and “exotic” molecules - those containing unusual (for pKa prediction) chemical elements such as transition metal atoms, and rare functional groups. There is

ACS Paragon Plus Environment


Page 2 of 24 2

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 course a separate and extensive research field of predicting protonation constants of side chains in proteins.32–35 In these systems a very large number of acidic and basic functional groups co-exist and this creates an entirely different challenge in comparison with all the challenges that we have mentioned above. In this work we will consider only molecules that have a small number (typically fewer than 10) of protonatable and deprotonatable functional groups. From the above discussion, it becomes clear that the number and the scope of issues facing a computational pKa predictor aspiring to be widely or universally applicable, are considerable. In the course of many years, we have been pursuing the development of a computational pKa predictor meant to be applicable to numerous types of practically important chemical structures. Our approach is based on the density functional theory (DFT)36 and works with three-dimensional representations of molecules. It is a complex workflow incorporating developments from the fields of cheminformatics, molecular mechanics, implicit solvation models, and, obviously, quantum mechanics. A set of empirical parameters used to partly account for close range effects that are not covered adequately by our solvation model and by our thermodynamic cycle adjust the so-called “raw” pKa values to those directly comparable to the experiment. A scheme that we call the “shell model” adjusts the use of these empirical parameters for different functional groups, and actually covers the whole chemical space, albeit with progressively decreasing accuracy for rarer and more exotic species. In this work, we review and discuss the latest developments of our program, the most significant of which are the shell model and an explicit accounting for multiple conformational forms. These improvements enable a much greater coverage of the chemical space and a higher accuracy and stability of pKa predictions. The progress is illustrated by a large number of diverse applications, empasizing drug-like molecules because those tend to be especially challenging owing to their larger atom count and greater conformational flexibility, and because an accurate pKa prediction of such molecules is known to be particularly valuable. We show that in certain cases the accuracy achieved can be so high as to correctly describe fine steric effects and correctly predict small differences in pKa of stereoisomers. The literature generally converges to the opinion that the accuracy targeted by pKa predictions should be below 1.0 pKa , and preferably as low as 0.5 units of mean unsigned error.14,16,37 Is this a reasonable aspiration? While it is claimed that accuracy for modern experimental pKa measurements lies within 0.1 pKa units,20,38,39 such assertions are made when referring to measurements conducted repeatedly in the same laboratory. When analyzing experimental pKa data published in the literature, we often observe that experimental measurements conducted on the same compound in different laboratories produce discrepancies as large as 0.5 pKa units, even in recent publications (see, for example, supplementary information in Ref. 39). We believe that differences in measurement techniques, sample purity and sample preparation protocols, solvent purity and content, and even temperature may all contribute to the observed variations. Because the un-

derlying parameterization and comparison of the final results is performed with experimental pKa values collected from various experimental groups, it would be unrealistic to expect to consistently achieve an accuracy of less than about 0.5 units with theoretical pKa predictions, regardless of the inherent accuracy of the theoretical methods or protocols used. This paper has the following organization. The section that comes next provides the basics of the DFT-based pKa prediction method. In Sections III-IV we discuss the main recent computational innovations of our program – the shell model and the treatment of conformations, respectively. Then, in Section VI, we show applications and benchmarks, and wrap up the discussion with a summary in Section VII.

II. METHOD

The basic methodology and the computational details underlying our pKa prediction method have not changed since the work in which the method was introduced.16 Here we briefly review the main points of that methodology and concentrate on the improvements that have built upon the original implementation over the years. All the scientific and technological improvements described below have been implemented in the computational program Jaguar pKa40 which uses the program Jaguar40,41 as its DFT engine. In the present section we assume that there is only one conformation in each protonation state. Our method is grounded in the thermodynamic cycle42 shown in Figure 1, which permits expression of the soughtafter free energy of acid dissociation ∆Ga through other free energies which can be either computed (∆Gg , ∆GAH sol , H+ A− ∆Gsol ) or taken from experiment (∆Gsol ): +

A H ∆Ga = ∆Gg − ∆GAH sol + ∆Gsol + ∆Gsol . −

(2.1)

The meaning of the free energies from this equation is made clear by Figure 1. The quantities pKa and ∆Ga are linked by the fundamental thermodynamic relation: pKa =

∆Ga . 2.303RT

(2.2)

This thermodynamic cycle has been used in numerous other works dedicated to pKa predictions (see, for example, publications 2,37,43 and references therein). Since quantum chemical methods are not designed to compute properties of sys+ tems without electrons, the experimental GH sol value of – 42 259.5 kcal/mol is utilized. A number of different and allegedly more accurate values for the solvation of the proton have been reported, (see, for example, Refs. 44,45 and references therein), and there might be a reason for opting for one of these values instead. However, choosing a different + value of GH sol would not change the results produced by our workflow because in our workflow the computed ∆Ga , which + depends linearly on ∆GH sol , would be effectively adjusted by linear empirical parameters for the generation of the final results (vide infra), and thus would merely translate into an ad-


Page 3 of 24

Journal of Chemical Theory and Computation 3

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

AH(g)

ΔGg

ΔGAHsol

AH(aq)

A–(g) + H+(g) ΔGA

ΔGa

– sol

ΔGH

+ sol

A–(aq) + H+(aq)

FIG. 1: The thermodynamic cycle used in our Jaguar pKa program. The subscript (g) denotes species in the gas phase while the subscript (aq) denotes species in the aqueous phase.

justment of these empirical parameters. Therefore, we keep + to the value of GH sol employed in Ref. 16 for consistency A− with this earlier work. The values ∆Gg , ∆GAH sol , ∆Gsol are evaluated in our protocol using DFT methods. The latter two free energies are evaluated with an implicit solvent model which solves the Poisson-Boltzmann equation.46–48 Experience shows48 that pKa values coming directly from equations 2.1, 2.2 can be a few pKa units away from the experimentally measured values, even for the simplest, most rigid molecules. Inherent deficiencies of DFT, basis sets, and the implicit solvation model are blamed for these large errors, which unfortunately cannot be corrected for all types of functional groups simultaneously by replacing or modifying these underlying components. For instance, when we replaced the B3LYP functional, standard to our pKa computation protocol, by the M06-2X functional49 (which has been shown to yield more accurate energetics in many cases)49–51 and reoptimized all the dependent parameters, this has not produced any systematic improvement to the root mean square deviations (RMSD’s) of the parameter fitting and to final pKa values across a large spectrum of functional groups (data shown in Supplementary Information). A good way of convincing oneself of the impracticality of using equation 2.2 directly for all functional groups is to realize that at 298K a mere 1.36 kcal/mol of error in ∆Ga translates into 1 pKa unit of error on the logarithmic scale.1,2,40 It is not to be expected that one can achieve an error significantly lower than 1.36 kcal/mol in calculation of ∆Ga , using DFT methods, across all types of functional groups and structures, no matter what thermodynamic cycle is used.2 This is because errors of bond dissociation energies produced by DFT methods and solvation energy errors are often above one and sometimes over several kcal/mol.45,52–55 Thus, equation 2.2, used directly, might yield accurate pKa predictions (with an error as low as 0.1 pKa unit) due to a a cancelation of errors or a special choice of the + value of ∆GH sol . However, this typically only works for one or a handful of functional groups at a time, with all other parameters of the computational protocol being equal.14,15 Explicit water molecules have been used in an attempt to improve the accuracy of pKa prediction for simple molecules,56,57 but it remains to be worked out how to apply such protocols for all types of chemistries in an automated fashion. We would like to emphasize that, in principle, a high accuracy in pKa prediction might be also achievable with a fully non-empirical ap-

proach, if one uses some of the largest basis sets, the most accurate levels of theory and solvation models (possibly in combination with adding explicit water molecules), together with a proper entropic treatment; however, such approaches would be impractical for all but the smallest and simplest chemical structures, due to their very high computational cost. However, recent works that use novel or sophisticated treatments to achieve high accuracy are worth paying attention to.33,58–65 The errors produced in computation of bond dissociation energies and solvation energies in our protocol have been shown to be systematic and to depend linearly on the functional group.16 This observation permitted development of a practical scheme that corrects the “raw” pKa obtained directly from equation 2.2 using a pair of linear parameters a, b optimized for each functional group: pKafinal = a pKa + b.

(2.3)

This optimization requires a training set consisting of sample molecules with experimental pKa values for each functional group. As a result, the authors of Ref. 16 achieved a dramatic reduction of pKa errors across a large area of the practically important chemical space. The idea of adjusting preliminary pKa values through a linear regression has been used in other pKa prediction protocols.66–68 The workflow’s quantum chemical calculations are performed as follows. Geometry optimizations of the protonated and deprotonated forms are carried out in either the gas or the solution phase (following user’s setting) with the B3LYP/631G* level of theory. Once the optimal geometries are available, accurate single point calculations are performed in gas and solution phases using the B3LYP functional, with ccpVTZ+ on the deprotonated atoms which bear formal negative charge and cc-pVTZ on the rest of the atoms. In evaluating solvation energies, our program uses the PBF implicit solvation model46–48 which has been shown to produce accurate solvation energies of neutral and ionic species. This solvation model solves the Poisson-Boltzmann equation explicitly, using a set of specially designed atomic radii that are vital for generation of quality molecular surfaces. The final solvation energies are corrected by a set of empirical parameters based on SMARTS patterns (see Section III for a discussion of the shell model, which serves as a mechanism for selecting appropriate empirical parameters). These parameters compensate for the first-shell interaction with the solvent which is missing



Page 4 of 24 4

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

FIG. 2: An illustration of the macro-pKa , micro-pKa , and nano-pKa terminology. Macro-pKa describes and encompasses (possibly) numerous protonation channels. Micro-pKa describes a single protonation channel characterized by (possibly) numerous conformations. Nano-pKa describes the protonation equilibrium between a single protonated conformation and a single deprotonated conformation. When micro-pKa ’s and macro-pKa ’s take multiple conformations into consideration (and are thus defined through nano-pKa ’s), they are referred to as conformationally-aware micro-pKa and macro-pKa , respectively.

from the implicit solvent model, and are absolutely necessary for generating accurate solvation energies for ionic species. The single point energies go into constructing the values of A− ∆Gg , GAH sol , Gsol via the following approximate formulas: − 5 ∆Gg ≈ EgA − EgAH + RT − T ∆S(H + ), 2

(2.4)

AH ∆GAH sol ≈ Esol ,

(2.5)

A ∆GA sol ≈ Esol .

(2.6)

−

−

In formula 2.4, the terms (5/2)RT and T ∆S(H + ) are the enthalpic and entropic contributions for H+ , respectively, under ideal gas conditions. The first contribution is equal to 1.48 kcal/mol, and the second is 7.76 kcal/mol, assuming room temperature and 1 atm of pressure. These are the standard approximations and values, and they have been used in

a number of other works.14,37,69,70 In all three equations 2.42.6 we neglect the entropic contributions for the species AH and A− , assuming that these contributions cancel out or at least can be partly absorbed by the empirical corrections from formula 2.3. This assumption, deemed reasonable for conformationally rigid molecules, is reconsidered in other context in Section IV. We also avoid computing zero point vibrational energies (ZPVE), and therefore enthalpies, in all our calculations because ZPVE are computationally expensive to generate. The overall ZPVE contributions to ∆Ga do not cancel out, but they are approximately constant for the same functional group16 and are therefore conveniently included in the empirical parameters a, b. Within a parameterization for a single functional group, the parameter a can be interpreted as compensating for the deficiencies in the variable components of equation 2.3, namely solvation and gas phase energies, whereas the parameter b effectively absorbs any innacu+ racies in the postulated constants GH sol = −259.5 kcal/mol and T ∆S(H + ) = 7.76 kcal/mol, as well as the approxi-


Page 5 of 24


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

mately constant vibrational and entropic terms neglected in equations 2.4-2.6. Now that all the free energies are computed, the final pKa can be assembled using equations 2.2-2.3 and then output for each set of the matched parameters a, b. There are approximately 100 specific functional groups (the exact number depends on the way one divides some of the groups into subgroups) in the current version of the Jaguar pKa program for which independent parameters a, b have been trained. Figure 3 shows these functional groups. The total number of experimental pKa data points encompassing all these groups is 1134. Note that in some cases functional groups “overlap”; that is, either the group is a subgroup of another group, or the compound can be thought of as belonging to both groups without clear precedence. As an example of the first case, consider a substituted 4-hydroxycoumarin which belongs to both the 4-hydroxycoumarin group and to the Het-OH group (a heterocycle linked to hydroxyl), with the first being more specific and therefore preferable. As an example of the second case consider 8-quinolinethiol which will match both the Ph-SH and Het-SH patterns, with it not being clear, at least without inspection of the respective training sets, which pattern one would prefer. In practice one can use an algorithm to choose the “best” set of parameters for a given compound. In the following section we will describe one such algorithm actually used in Jaguar pKa, although alternative algorithms should not be hard to formulate. The discussion in the above two paragraphs makes it clear that an association of the predicted pKa with the protonation/deprotonation of a specific functional group in the given molecule plays a major role in our approach with such a value generated for a single process thus being a micro-pKa (for the distinction between micro- and macro-pKa see Figure 2). If the molecule does not contain other protonatable/deprotonatable groups, the micro-pKa predicted by the program should correspond directly to the macro-pKa observed experimentally. If several protonatable/deprotonatable functional groups are present, their micro-pKa ’s can be computed independently and then combined into a macro-pKa using an elementary algebraic equation. In practical calculations often only one functional group should be considered, even if several protonatable/deprotonatable groups are present in the molecule. This is because it can be obvious from the chemical intuition which functional group is primarily responsible for the observed macro-pKa . For instance, if the molecule contains both phenolic hydroxyl and aliphatic carboxyl groups, and we are interested in predicting the pKa1 constant, it is clear that we should only compute the micro-pKa of the carboxyl group. The micro-pKa of the phenolic hydroxyl group should be several logarithmic units away from that of the aliphatic carboxyl group, and as such provides a totally negligible contribution to the total effective pKa1 . If there are two or more symmetrically-equivalent protonation/deprotonation sites in the molecule, there is no need to compute the micro-pKa of all of these sites for their subsequent inclusion into the formula for macro-pKa . Computing micro-pKa for just one such group will suffice, upon the conjecture that the micro-pKa of the rest of the groups is the same. Jaguar pKa can identify symmetrically-equivalent

sites automatically by analyzing the molecule’s SMARTS pattern,71 compute the micro-pKa of one of them, and report the macro-pKa , accounting for the micro-pKa ’s of the symmetrically-equivalent groups that it did not compute explicitly. This saves computational time and frees the user from the burden of indentifying symmetry “by eye” and assembling the macro-pKa “manually”. In the existing implementation of Jaguar pKa, the user must decide which protonation/deprotonation process is being considered in the pKa prediction and then pick the functional groups thought to be involved in this process. This decision making process should be straightforward to automate in a future version of Jaguar pKa. When the functional groups are selected the program automatically identifies them, matches them against the SMARTS patterns representing the functional groups from Figure 3, and in this way selects the empirical parameters a, b for use at the end of the calculation. Then it generates the corresponding protonated and deprotonated species and launches all necessary calculations. Conformational treatment of the generated species is optional, although greatly recommended, and is discussed in Section IV.

III. SHELL MODEL

As we have seen in the previous section, an application of empirical parameters via Equation 2.3 is necessary for accurate prediction of pKa . The selection of parameters a, b, specific for each functional group, is carried out with the help of the shell model which is elucidated in this section. The shell model is a tree-like graph that organizes SMARTS patterns, and possibly different types of data associated with each of these patterns, in a hierarchical order of “shells”, which are essentially the graph’s nodes. The shell model does not have to be a tree, but in application to our pKa prediction scheme it usually is, and we will assume it to be such in the rest of the text. The tree is organized from top to bottom, with the SMARTS patterns ranging from most generic (at the top) to very specific (at the bottom). The presence of the most generic possible pattern guarantees that any given molecule will match at least one pattern. With the hierarchical organization in place, it is possible to automate the classification of a given molecular structure in terms of its association with a specific type of chemistry. Once the association is established, one can automatically apply settings and parameters belonging to the matched pattern. Therefore, one passes a structure to the computer program as an input and obtains, as an output, settings and parameters tailored for the treatment of this particular kind of molecule. An shortened version of the shell model, together with its application to the classification of a benzodiazepine protonation, is shown in Figure 4. The actual shell model used in the current version of Jaguar pKa is more extensive, as it comprises over a hundred shells based on the patterns from Figure 3, and would be difficult to put on a single page. It is, however, available in the Supplementary Information. If the model from Figure 4 is used in conjunction with the pKa prediction described in Section II, then each shell shown in the



Page 6 of 24 6

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

O

O

O

O

N H

N H

H

N

O

N

H

H N+

O-

H

N

N

N

OH

OH

H N+

OH

H

N

N

N

N H N R

N H

R'

H N

N

N H R

R

R N H

C

NHOH

OH Het

N H

NRR'

C

O

R

Ar OH

R

O

S

N Het

Het Het

OH

R R R'

O

O R

R

S

OH O

N

NH O

N

O R' S R'

R O

N

OH R

Alk OH Alk

OH

N

O

S

N

R O

R

Alk

R

Alk

N

SO2R'' R''

OH R

R'' N H

R' N

N R'RC

R''

N H

B OH SH

OH

N H NH2 NHR SH

R''

R'

OH

H C

H N

P

Ar

OH

N

OH

O

P

R'

Alk

N

NRR'

R'RN

NR''R'''

N

N

N

OH

RR'R''N RR'R''P RR'S

N H

N H

+

N

N

R

OR'

OH

O

O OH

N

O-

R''

OR OH

R'

+

O O

R''

H H

NR'''R''''

N H

OH

P

H

R

O

RO P

O

CONRR'

Ar

R

OH

OR'

COR

Se

N

R'

O

RO P

O

O

O

O

R

N H

R

R'

R' H N

R

R'

Het C

R

CH

O

Ar

N

CH

NR'R'' R''

NR'''R''''

R' N

R''' O

R

H

CN

OH

OH

R''

O

R''

OH R

O O

R

R

R'

N H

O

H

O

O

R

O- R As OH R As OH R' OH

Ar +N

N OH

R'

N H

OH

O

Ar

R

N H

NRR'

O

N

N H

O

O

N

N NH

NH

+ G O

O- RR'R''N RR'R''P RR'S

G

Ar

O-

H H N N Ar N N Alk O O H R Ar Ar' C S R S R' O S OH S S O S SH R R O O O R R' H H H H R'' O N N N N R R NR'R'' NR'R'' N R R' Alk R OOH R R N R'' S S S O O N O O H H OH H NH O N O N R NO2 S OH OH HN N R H O R' R' NO2, NO OH

O

H N

FIG. 3: The functional groups that have an explicit parameterization in Jaguar pKa. The hydrogen atoms marked in red indicate the sites of deprotonation, while the atoms marked in blue indicate the sites of protonation. Alk denotes alkyl groups, Ar denotes aryl groups, and Het stands for heterocycles. G is any group of atoms and R typically begins with a carbon atom. Unmarked H atoms can be replaced by essentially any substituents without affecting the recognition of the functional group.

image has a SMARTS pattern and parameters a, b associated with it. The shells marked in gray are those that matched the benzodiazepine functional group. In our implementation, the shells’ SMARTS patterns correspond to protonated functional groups and are matched to the structure that has already been protonated. We chose the protonated patterns as the greater number of atoms leaves less room for ambiguity when matching. Each set of parameters a, b corresponding to the matched shell leads to its own pKa prediction using Equation 2.3. Of

course, as the matching becomes more and more specific, as when going from the top generic atom shell to the benzodiazepenes shell in Figure 4 the parameters a, b are deemed to be progressively more relevant, and therefore more accurate. So it is probably the pKa prediction corresponding to the most specific pattern that the user is going to prefer in most cases. Now it is appropriate to discuss the parameterization of the shells. All shells can be classified into elementary and derivative. An elementary shell is the shell that is an end vertex of


Page 7 of 24


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

FIG. 4: The shell model that was used for classification of a protonation of benzodiazepene derivative. Only a part of the model could be displayed in this image. The shells that could not be displayed, for the economy of space, are marked by the ”...” symbol. The full model is available in the Supplementary Information. The nitrogen atom to be protonated in the benzodiazepene molecule is shown by the gray circle. The shells marked in gray are those that match the protonated benzodiazepene functional group.

the tree-like structure; it is directly connected only to shells that have a less specific SMARTS pattern. All shells that are not elementary are derivative. Thus, in Figure 4 “sulfonic acid”, “benzodiazepine”, and “conjugated carboxylic acid” are examples of elementary shells and “sp3-like O”, “imidazole as acid”, and “imine” are examples of derivative shells. Only elementary shells need to be explicitly parameterized through specific training set data points in our implementation. Derivative shells obtain their parameters a, b through a simple automatic model that extends their fitting to all of the training set points for the underlying elementary shells. Again, using Figure 4 as an example, parameters of the “imine” shell are constructed from the training set data for the elementary shells “benzodiazepine”, “hydrazone”, and the other elementary shells belonging to “imine”, but not explicitly shown in the Figure. The parameters of the higherlying “sp2-like aliphatic N” shell are constructed from all the training set data of its underlying elementary shells, namely “barbituric acid”, “hydrazone”, “benzodiazepine”, and other elementary shells grouped under “sp2-like aliphatic N”, and not explicitly shown in the Figure. The parameterization of an elementary shell is carried out as follows. One chooses several chemical structures (typically more than four and fewer than twenty) that satisfy the following criteria: (i) they contain a single functional group matching the elementary shell’s SMARTS pattern (ii) their experimental pKa ’s are known with a high degree of trustworthiness, measured in water, at or around 298K, and preferably in absence of salts and other solvents (iii) the structures contain either a single functional group or several functional groups, but in the latter case under

conditions where the other groups are certainly not interfering with the protonation/deprotonation process of the principal functional group (iv) the structures are preferably simple and non-flexible (v) the structures undergo no tautomeric equilibria in their solutions If several similar experimental pKa values exist for the given structure they are averaged, and the result is associated with the structure. The idea behind selection rules (i)-(v) is that we would like to attribute the measured pKa of the compound to the protonation/deprotonation processes around one functional group, free of other complicating factors such as effects from other functional groups or the presence of alternative conformational and tautomeric forms. Once the compounds for the ‘training’ set are prepared in accordance with rules (i)-(v), their “raw” pKa ’s are computed following the thermodynamic cycle from Figure 1. The parameters a, b are the coefficients of the linear fit of the raw pKa values to the corresponding experimental values. Typically, there is good to very good correlation (the r2 coefficient is 0.8 and higher), at which point the parameters a, b are associated with the shell and the parameterization is declared complete. However, there are sometimes cases when the correlation is poor. This is a sign that either there is something wrong in the interpretation of the experimental results or in the experimental results themselves, or that the physics of the protonation/deprotonation process is more complex than we have heretofore assumed. In such cases we carefully re-examine the experimental pKa ’s. Then, if we still have no reason to doubt the correctness of the experimental measurements, we split the shell into two or more other sub-shells, possibly adding more data points, and repeat the parameterization procedure for each of the sub-shells, until the satisfactory linear correlation is obtained for each of



Page 8 of 24 8

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 sub-shells. The shell which is a “parent” of the sub-shells thus becomes a derivative shell and its parameters a, b are obtained simply from a linear fit of all the individual data points comprising its elementary shells. The correlations thus obtained for derivative shells are necessarily worse than correlations obtained for elementary shells. This is why pKa predictions produced by the parameters of the matched elementary shells are preferable to those generated by the parameters of the matched derivative shells. Derivative shells are, however, indispensable for coverage of all chemical space because there may be cases when a given compound is not matched by any of the elementary shells present in the current shell model. However, a match might be achieved by one of the derivative shells, and certainly by the most generic shell, which matches all compounds. As in the situation above, the pKa value coming from the most specific derivative shell must be preferred. The parameterization procedure is somewhat laborious and not necessarily objective. The person who performs the parameterization must decide whether the training set should include a particular data point or not, and also how the shell will be split into sub-shells. In the future, we would like to work toward an implementation of a data-scouting computer algorithm that would ease the burden of the parameterization on the developer and perform the procedure in a more objective way.

IV. TREATMENT OF CONFORMATIONS

As the original paper on Jaguar pKa indicated,16 the conformational flexibility of molecules is a major issue for quantum chemical pKa prediction. This is because the micro-pKa ’s of individual conformations might be dispersed across a broad range of values, and preferring any particular conformation over the other conformations without taking the corresponding energies into account is unjustifiable. The results shown in the following section clearly demonstrate that the accuracy of pKa prediction in the absence of any special treatment of the conformational problem depends strongly on the conformational flexibility of the molecule and on its input conformation. In this section we discuss how this problem can be addressed, and describe a solution recently implemented in Jaguar pKa. The simplest way to deal with multiple conformations is to assume that only one conformation, identical for the protonated and the deprotonated forms, is dominant, and is therefore responsible for the observed pKa . The next level of approximation is to recognize that the protonated and the deprotonated forms of the same species might have different dominant conformations, as is expected for structures that are capable of forming intramolecular hydrogen bonds. Ultimately, one can accept that more than one protonated and more than one deprotonated conformation might exist in appreciable quantities in solution, and all these species can contribute to the observed pKa . Whichever of these approximations is adopted for practical realization, there should be an independent mechanism to select one or more significant conformations prior to subjecting them to quantum chemical

calculations. Since there are numerous conformational sampling methods with settings potentially resulting in different sets of conformations, one can in principle have a large number of ways to implement conformationally-aware quantum mechanical pKa predictions. Early versions of Jaguar pKa approached the conformational problem from the conceptually and computationally more accessible side: by assuming that only one conformation per protonation state contributed to pKa . The most energetically favorable protonated and deprotonated states were selected by the associated molecular mechanics program MacroModel.72 This program employs the OPLS2005 force field73–75 with settings that ensure a reasonable compromise between the breadth of the CS and its computational cost. A new and more accurate version of the force field, OPLS3,76 has not been yet used in Jaguar pKa systematically. An investigation into the possible improvements for pKa prediction results due to a transition from OPLS2005 to OPLS3 is reserved for a future work. Our numerous applications of this methodology to large and flexible molecules convinced us that in some cases such a basic treatment of conformations was inadequate for achieving an acceptable error of less than 1 pKa unit. Another problem that we observed was that the final pKa prediction depended, in some cases significantly, on the initial choice of the conformation that was passed to MacroModel. It was possible to target a higher accuracy by generating conformations with a custom protocol, prior to performing the pKa prediction itself. One could then even statistically average the pKa ’s computed for individual conformations. However, such practices were not standardized, and their effects on the final results have not been systematically investigated. Clearly, an improvement to the handling of conformations was needed, and a more thorough study of the influence of the CS settings on the accuracy of pKa prediction was desirable.

A.

Traditional macro-pKa

In order to derive and explain the formulas that we use to compute conformationally-aware macro-pKa values, it is best to start with the well-known formulas for computing macropKa from micro-pKa ’s, assuming that each proton dissociation channel is represented by a single protonated and a single deprotonated form. Then we generalize this formula to multiple conformations, which enables us to compute a macro-pKa involving multiple protonation channels with any number of conformations. Suppose a molecule in its aqueous solution exists as a set of “protonated” tautomers {P} composed of N species P1 , P2 , ..., PN , each bearing total charge C, and a set of “deprotonated” tautomers {D} composed of M species D1 , D2 , ..., DM , each bearing total charge C − 1. For the simplicity of discussion we will assume, without losing generality, that the ordering of these two sets of protonation states goes from the most populated to the least populated forms. Other protonation states with charges C +1, C −2, etc., might exist but we are not concerned about them here because our


Page 9 of 24


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

goal is to compute the macro-pKa for the transition from C to C − 1. Once this formula is found it can be equally applied to other charge transitions. By definition, for the dissociation equilibrium P ⇋ D + H+

(4.1)

we have Ka =

[D][H+ ] ([D1 ] + [D2 ] + ... + [DM ]) [H+ ] = , (4.2) [P] [P1 ] + [P2 ] + ... + [PN ]

where [P] and [D] are the total equilibrium concentrations of the forms for P and D, respectively. Expanding this expression into M terms with the numerators [D1 ][H+ ], [D2 ][H+ ], and so on, and then dividing both the numerator and the denominator in each of these terms by the numerator, we cast equation 4.2 into an actionable form Ka =

M X j=1

PN

1

1 i=1 Kij

,

of conformations. However, formula 4.3 ceases to be valid for conformations as soon as we replace Kij in it by the corresponding micro-pKa ’s computed on the sets {P} and {D} with the thermodynamic cycle from Figure 1 (or, in fact, with any other similar thermodynamic cycle or empirical formula used in practice). To convince oneself of the inapplicability of the naive interpretation of formula 4.3 for conformations, one can consider a simple case of the dissociation of phenol. The deprotonated form is represented by only one conformation (M = 1), whereas the protonated form is represented in principle by an infinite number of conformations characterized by the differerent dihedral angles CCOH. There is a continuum of these protonated conformations, all having similar energy. Note that nothing compels us to consider only conformations occupying a minimum on the potential energy surface. Since all the conformations have similar energy, their Ki1 constants will be similar, leading to a nonsensical result in which Ka tends to zero (and the respective pKa to infinity) with the increasing number of conformations considered.

(4.3) B. Conformationally-aware macro-pKa

where Kij =

[Dj ][H+ ] [Pi ]

(4.4)

is the micro-dissociation constant for the deprotonation channel Pi ⇋ Dj + H+

(4.5)

Micro equilibrium constants 4.4 can be readily computed by Jaguar pKa. Importantly, equation 4.3 “converges” with the number of dissociation channels. It is easy to see from equation 4.3 that negligible dissociation channels, characterized either by very low concentrations of the protonated form or the deprotonated form, will contribute negligibly to Ka . So that for practical purposes we can effectively truncate the summations in equation 4.3 by introducing the “effective” number of deprotonated and protonated forms, Me and Ne , respectively. Obviously, failure to include highly populated forms from the sets Pi or Dj , or both, in the summations can lead to large errors in prediction of Ka . Equation 4.3 can be a working formula for empirical approaches that do not have to deal with conformations.20,77 It can also be used in quantum mechanical pKa calculations if we assume that the protonation states {P} and {D} are each represented by a single conformation.23 It might be tempting to use equation 4.3 for conformations in the same way one uses it for protonation states, by interpreting {P} as a series of protonated conformations and {D} as a series of deprotonated conformations. However, it is easy to see that the formula essentially breaks down and becomes disfunctional when applied to conformations. Equation 4.2 can still be applied to conformations, provided we assume that we can come up with a reasonable definition of a “concentration” of a specific conformation when there is a continuum

An appropriate way to treat conformations with our thermodynamic cycle would be to properly take into consideration the conformational flexibility of the molecules in each protonation channel by including the corresponding entropic A− H+ contributions in the terms ∆Gg , ∆GAH sol , ∆Gsol , and ∆Gsol in equation 2.1. This path would require fundamental and non-trivial changes in the construction of our workflow, and therefore we decided to pursue an alternative and conceptually much simpler approach that would still allow us to treat the problem of the conformations. Let us assume that we are dealing with a single protonation channel, that is, there is only one protonated tautomer and one deprotonated tautomer, each characterized by a number of conformations. Once we discuss the situation with a single protonation channel we will extend the conclusions to multiple protonation channels. For the clarity of the following discussion it appears apposite to introduce the notion of the nano equilibrium constant (and the corresponding nano-pKa ’s), which is the equilibrium constant characterizing a protonation equilibrium between a single protonated conformation and a single deprotonated conformation. Thus, nano-constants characterize protonation equilibria between individual conformations in a single protonation channel. Micro-constants describe a single protonation channel (possibly) involving multiple conformations. Finally, macro-constants describe (possibly) multiple protonation channels and should be directly comparable to equilibrium constants derived from the experiment. The pKa from formula 2.3, computed using the thermodynamic cycle that involves a single deprotonated conformation and a single protonated conformation, and further adjusted by the empirical parameters a, b, could serve as a definition of nano-pKa . However, as we have seen, the nano-pKa defined through formula 2.3 would have a serious drawback of not accounting for the entropy of the conformations, es-



Page 10 of 24 10

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

pecially in treatments that consider conformations explicitly. The missing entropic prefactor e∆Sconf /R can be actually included through the following Ansatz: nano Kij =

λDj final K . λPi ij

(4.6)

final Here, Kij is the equilibrium constant extracted from final pKa , which is defined in Equation 2.3. Next, λDj is the probability of finding the deprotonated conformation Dj in the conformational ensemble {D} comprising M deprotonated conformations. Similarly, λPi is the probability of finding the protonated conformation Pi in the ensemble {P} comprising N protonated conformations. The λ factors account for the different populations of the deprotonated and the protonated conformations associated with the nano-constant. We will use Equation 4.6 as a definition of the nano equilibrium constant. The corresponding nano-pKa is obviously computed by taknano ing the negative logarithm of Kij .

One way to define the probabilities λDj and λPi is through the conformational energies: λDj = e−EDj /kT /

M X

e−EDl /kT ,

(4.7)

The summation limits M conf and N conf are the number of deprotonated and protonated conformations, respectively. Formulas 4.11 and 4.12 are particularly convenient in practice as λDj and λPi thus defined are computed directly from pKa values output by the Jaguar pKa program, with the empirical corrections included. We tested both definitions of the λ factors in some of our calculations and concluded that the final pKa ’s produced by both approaches are very similar. We chose equations 4.11, 4.12 as our working formulas, and therefore they were used in all the multiconformational pKa calculations reported in this work. Now that formulas 4.6, 4.11, and 4.12 can be used for computing the nano equilibrium constant for a single protonation channel, we can move on to define the conformationallyaware micro equilibrium constant for a single deprotonation channel, which should be computed from nano equilibrium constants. For this purpose one can adopt the methodology expressed in relations 4.1-4.3, this time viewing P1 , P2 , ..., PN not as protonated tautomers but as protonated conformations, and, analogously, viewing D1 , D2 , ..., DM not as deprotonated tautomers but as deprotonated conformations. This permits us to define the conformationally-aware micro-dissociation constant:

l

λPi = e

−EPi /kT

/

N X

Kamicro e

−EPl /kT

.

Here, k is the Boltzmann constant, T is the absolute temperature, and EDj and EPi are the energies of the conformations Dj and Pi , respectively, computed in either the gas or the solution phase, depending upon the particular computational protocol. An alternative way starts with concentrations: λ Dj =

[Dj ] , [D1 ] + [D2 ] + ... + [DM ]

(4.9)

λPi =

[Pi ] . [P1 ] + [P2 ] + ... + [PN ]

(4.10)

The numerators and the denominators of formulas 4.9 and 4.10 can be multiplied by [H+ ]/[P1 ] and [D1 ]/[H+ ], respectively, and after some algebraic rearrangement and equating pKa with pKafinal , the expressions for λDj and λPi become directly connected with the empirical corrections from formula 2.3: λDj = 1/

conf M X

final

10pKa1j

−pKafinal 1l

,

(4.11)

.

(4.12)

l

λPi = 1/

conf NX

l

final

10pKal1

−pKafinal i1

j=1

(4.8)

l

=

conf M X

1 PN conf i=1

1 nano Kij

.

(4.13)

Here, M conf and N conf , as above, refer to the number of deprotonated and protonated conformations, and the summations go over conformations. Notice that owing to the λ factors the generalized formula 4.13 behaves correctly regardless of the number of conformations one takes. In the case of a single protonated and a single deprotonated conformation in each protonation channel, 4.13 is equivalent to the earlier formula 4.3. When there are several conformations considered for either protonated and/or deprotonated forms, the nano equilibrium connano stants Kij are properly weighed by the entropic λ factors. This prevents the kind of collapse that formula 4.3 suffered from, when used earlier for conformations in lieu of tautomers. As a simple example of the correct behavior of the conformationally-aware formula 4.13, consider Q (possibly symmetrically equivalent) conformations with exactly the same energy making up the deprotonated form (M conf = Q), and a single protonated conformation (N conf = 1). In such a case the earlier formula 4.3 for Ka , used for conformations, reduces to K11 + K12 + ... + K1Q , where K1j are all the same because of the isoenergeticity of the deprotonated conformations. This result is nonsensical because Ka depends on the number of equivalent conformations Q. The generalized and conformationally-aware formula 4.13, however, does not have such problems, as Kamicro reduces to nano nano nano K11 + K12 + ... + K1Q . The λ-factor of each isonenergetic deprotonated conformation is λDj = 1/Q, irrespective of whether one uses formula 4.7 or formula 4.11, leading to nano final final K1j = (1/Q)K1j . And, as all constants K1j are equiv-


Page 11 of 24


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

nano nano nano final alent, the sum K11 + K12 + ... + K1Q turns into K11 , micro independent of the number of considered rendering Ka isoenergetic conformations Q. In practical calculations one needs to make a decision about which sets of conformations {P} and {D} to use. The Jaguar pKa program lets the user choose the settings for the MacroModel component, controlling the maximum number of conformations and the energy window from which these conformations will be selected. For the selection of more relevant conformations we also relaxed a number of constraints originally used in the MacroModel conformational search for reasons of efficiency. This increased the overall computational time slightly, but resulted in higher accuracy of the final pKa predictions, as the results in the following section indicate. Once the conformations are selected by MacroModel they are passed to Jaguar pKa for further refinement in DFT geometry optimizations. At this point the user has a choice of whether to optimize the initial structures coming from MacroModel in the gas phase, which is computationally less expensive, or in the solution phase. The following section compares the quality of results obtained with both of these types of optimization. Finally, it is time to consider multiple protonation channels and the corresponding formulas for the conformationallyaware macro-pKa . This task is trivial, because the formula for the macro equilibrium constant has already been presented in Equation 4.3, which now should use conformationally-aware micro for each protonation channel involvmicro-constants Kij ing dissociation of the protonated tautomer Pi into the deprotonated tautomer Dj :

Kamacro

=

taut M X

j=1

1 PN taut i=1

1 micro Kij

.

(4.14)

Here, M taut is the number of deprotonated tautomers and N taut is the number of protonated tautomers and the summations go over tautomers. Formula 4.14 in conjunction with the conformationally-aware micro-dissociation constant 4.13 is a generalization of the well-known formula 4.3. To summarize, this section defines formulas for computing the conformationally-aware macro-pKa from conformationally-aware micro-pKa ’s (formula 4.14). The conformationally-aware micro-pKa is defined through nano-pKa ’s (formula 4.13), and the nano-pKa ’s, in their turn, are defined through relation 4.6 and are computed from a thermodynamic cycle that involves a pair of individual conformations at a time, with further corrections compensating for the neglect of explicit solvation effects and entropic contributions to free energies.

V. PARALLELIZATION

Jaguar pKa calculations, especially those involving larger molecules and several conformations, can be computationally expensive. Multiple individual DFT calculations, including geometry optimizations, are performed in every pKa prediction. Geometry optimizations in the solution phase are par-

ticularly time-consuming. Therefore, it is important that one should be able to perform a Jaguar pKa calculation on multiple central processing units (CPU). Many individual DFT calculations in our computational protocols are completely independent (especially when acting on multiple conformations), which helps achieve good distributed parallelization. Conformational search, which precedes the DFT calculations, can also be performed in parallel using the standard MacroModel options. Recently, we have been able to complete distributed Jaguar pKa calculations maximizing the number of individual DFT calculations, with with up to 100 CPUs utilized in one pKa prediction. As an example of timings relevant to practical calculations consider Jaguar pKa calculations performed on the drug atenolol having molecular formula C14 H22 N2 O3 . A protocol involving 10 conformations, gas phase geometry optimizations, and an aggressive conformational search, applied to this molecule took approximately 1 day, 5 hours, and 2 minutes to complete on 1 CPU from the AMD Opteron Processor 6274 chip. A calculation with exactly the same settings, except carried out on 10 CPUs, took approximately 3 hours 45 minutes to complete, showing a 7.7 times speedup. Finally, an analogous calculations but with 20 CPUs took approximately 2 hours 17 minutes to complete, resulting in a 12.6 times speedup. Our Jaguar DFT benchmarks on geometry optimizations and vibrational frequency calculations performed with the above-mentioned AMD processors and with Intel Xeon E5-160 processors show that the latter is up to a factor of 4.0 faster than the former, although the scaling with the number of CPUs remains approximately the same. We have not measured the performance of the Xeon processors in application to our pKa prediction protocols, but we expect the improvements of the timings by at least a factor of 2 to 3, for the equivalent number of CPUs, in the comparison with the above mentioned timings on the AMD chips. All timings given in this section refer to wall clock time.

VI. ILLUSTRATIVE RESULTS

A truly universal pKa predictor should achieve a good accuracy on diverse sets of functional groups and molecule types. While we acknowledge that for some special types of structures the results currently produced by Jaguar pKa are less accurate than desirable (see the discussion in Section VII), we believe that the accuracy achieved across the chemical space representing the larger and conformationally flexible drug-like molecules is very good. This section shows Jaguar pKa predictions on several different functional groups in drug-like molecules as a function of the number of conformations taken, the thoroughness of the preliminary selection of conformations by MacroModel, and the type of the subsequent geometry optimizations performed on this conformations with the DFT method. Most calculations consider a single protonation channel. Exceptions are symmetricallyequivalent protonation channels, which are accounted for automatically, and cases where we combined micro equilibrium constants for two protonation channels into a macro-constant



Page 12 of 24 12

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

FIG. 5: Differences in pKa values obtained using two distinct conformations A and B as input for cyclic amidines from Ref. 78. In plot (a) the x-axis refers to the structures listed in Table I, and the y-axis shows the pKa differences for different protocols. The ideal result for every structure would be zero difference in pKa prediction and the corresponding point on the plot would lie on the horizontal orange line. All protocols involving conformational search included solution phase optimizations. The “no conformational search” protocol (red line) used optimizations in the gas phase. The conformations A and B for all the structures are schematically shown in diagram (b) and have RMSD around 0.90 Angstrom and larger.

guanidines by the most sophisticated model likely borders on the level of the noise and capabilities of DFT to describe the dissociation equilibria.

using Equation 4.14.

A.

Cyclic amidines

One class of drug-like molecules the pKa ’s of which have received close attention in recent publications is cyclic amidines.78–81 These are typically β-secretase (BACE) inhibitors, and their pKa ’s lie in a broad range of biologically important values. In Table I we present detailed Jaguar pKa results on a series of sixteen cyclic amidines from Ref. 78. It is evident from this table that the results systematically and dramatically improve upon transition from the simplest protocol, in which no CS was performed, to the sophisticated protocol which involved an aggressive CS, up to 10 conformations per protonation form (M conf , N conf = 10), and geometry optimization in implicit solvent with the PBF method. Both the mean unsigned errors (MUE’s) and the number of outliers go down. For stereoisomers 66 and 67 marked by a substantial difference in their experimental pKa values equal to 0.7 units, Jaguar pKa ’ was able to correctly reproduce the individual pKa , as well as the experimental difference. Another pair of stereoisomers, 89 and 109, has a small difference in experimental pKa ’s equal to 0.3 units, which approaches the limits of the theoretical pKa prediction. We therefore do not consider it a problem that Jaguar pKa was not able to distinguish this pair of stereoisomers correctly (predicting an inverse difference of 0.3 units with the most accurate protocol). Overall, the extremely accurate result provided for this test set of

Besides the diminishing errors, another criterion that could show improvement in the protocols is a lower sensitivity of the results to the input conformation. Ideally, pKa prediction should not depend on the initial conformation specified by the user, and CS should find the lowest energy conformations starting from any input conformation. In practice, however, the result might be sensitive to the CS parameters used, and in some cases fail to arrive at a virtually identical set of conformations. Figure 5 shows that indeed, as we go from the least sophisticated protocol (no CS) to the most sophisticated protocol (aggressive CS and 10 conformations per protonation form), the pKa prediction results become less dependent on the initial structure. The aggressive CS (part of the protocols for the lower two plots in Figure 5) leads to dramatic improvement in the quality of results. From this figure it appears as though the number of conformations taken does not have a pronounced effect on the sensitivity to the initial structure, provided that aggressive CS is performed. This is true for this current test set, and is probably attributed to the fact that in this case both the protocols that involve the aggressive CS are already extremely accurate and approach the level of the numerical noise that would be hard (and unnecessary) to transcend. However, the results obtained on other test sets (vide infra) clearly indicate that the greater number of conformations in combination with aggressive CS leads to an increased stability of the results vis-à-vis the choice of the input conformation.


Page 13 of 24


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

A more detailed investigation of the importance of the conformational treatment was undertaken for three cyclic amidines from the work of Rombouts et al.80 Three structures, representative of the BACE inhibitors discussed in that work, were essentially randomly picked for assessment of our pKa prediction protocols. For each of the three structures, Table II shows standard deviations computed for a set of 15 different pKa predictions that used different and distinct conformations as input (RMSDs computed for all pairs of the relevant conformations ranged from 0.92 to 4.37 Angstrom). In an ideal protocol, pKa predictions would not depend on the initial conformation. Then all 15 calculations here would produce exactly the same pKa and the corresponding standard deviation would be zero. In reality, not all the pKa values are the same. The greatest variability in the pKa predictions is observed with the least sophisticated Jaguar pKa protocol, that is, the protocol avoiding the CS, with the subsequent geometry optimization in the gas phase. In this protocol, standard deviations as large as 2.092 were recorded for the (2S-3R)-7a molecule, which is the result of the scatter from pKa 0.8 to 8.1, depending on the initial conformation, while the experimental pKa is 7.6. The next least sophisticated protocol, the one that avoids the CS but performs geometry optimization with the implicit solvation model, shows a marked reduction in the standard deviations for all three molecules. A futher reduction of the deviation is seen for the protocols involving the aggressive CS. These results are in qualitative agreement with the results shown in Figure 5 and in Table I for related cyclic amidines. Even though 1 conformation was sufficient to achieve what appears to be convergent results in Figure 5 and in Table I, the same protocol with 1 conformation is not satisfactory in the more detailed investigation reported in Table II. Take for instance the molecule (2R,3R)-7a, for which the largest standard deviation was recorded in a 1-conformation protocol (0.563). For this molecule, one outlier with pKa = 10.0 was detected (the experimental value is 7.8), even though the remaining forteen predictions in the set lie close to one another and to the experimental value. Another example is the set of 15 pKa predictions for the molecule (2S,3R)-7a, that used a 1-conformation protocol involving a gas-phase optimization. Its significantly large standard deviation of 0.423 follows from the scatter of pKa values between 6.9 and 8.3 (while the experimental value is 7.6) that is common in this set. These large variations disappear when one uses 5 conformations per protonation species, regardless of the subsequent geometry optimization protocol, and give standard deviations that are modest and acceptable. For the largest 5-conformation protocol standard deviation of 0.091 [(2S,3R)-7a, gas phase optimization], the individual pKa values ranged from 7.2 to 7.6, depending on the initial conformation. The stability achieved for 5 conformations cannot be said to noticeably improve upon further extension to 10 conformations. This is especially evident in the set characterized by the lowest standard deviation of 0.032 (2R-6a, 5 conformations, gas phase optimization), where a remarkable stability was achieved with the variation of the predicted pKa values for 15 individual calculations using different input conformations lying between 9.5 and 9.6 (with the experimental pKa

being 9.6).

B. Tertiary amines

The tertiary amino group is a ubiquitous structural motif in drug-like compounds, partly owing to the variability of its micro-pKa near the blood brain barrier. In Table III we present Jaguar pKa predictions using a test set of large and flexible tertiary amines from Ref. 82. Two micro-pKa calculations were conducted for the protonation of each of the nitrogen atoms of the piperidine ring (compounds with A5 substituent), and then manually combined using formula 4.14 to produce the final macro-pKa ’s. The present test set of tertiary amines possesses MUE’s that are somewhat larger than those seen for cyclic amidines. As in the previous tests, the quality of the results (as inferred from the value of the MUE and the number of outliers) improves when performing a more aggressive conformational search, and when utilizing a greater number of conformations. However, geometry optimizations in the solution phase produced larger MUE’s than those obtained with gase phase optimizations. The number of outliers for geometry optimizations in the solution phase was also quite large and did not go down with an increase in the number of conformations taken, in contrast to geometry optimizations in the gas phase, where the trend was as expected. In all other tests of the conformational protocols that we performed (those shown and not shown in this work), including those on tertiary amines, geometry optimizations in solution produced results that were generally better or at least no worse than gas phase optimizations. Therefore we conclude that the somewhat disappointing performance of solution phase optimizations in Table III is an uncommon case. A very interesting class of tertiary amines from the point of view of their pKa properties, is the so-called proton sponges.83–85 These compounds display unusually high basicities as a consequence of a sterical strain leading to a repulsion of unpaired electrons on tertiary nitrogen atoms, which can be relieved by binding a proton. A proton sponge may possess a pKa of around 12, whereas a related sterically nonhindered amine may have a dramatically different pKa of around 5. DFT-based pKa predictors should excel for these types of compounds, in contrast to empirical pKa predictors which, failing to properly estimate or even recognize sterical strain, should see essentially no difference in pKa for sterically strained versus sterically non-strained compounds. Proton sponges and their pKa ’s have been studied with DFTbased methods before.86–88 Here, we have applied Jaguar pKa to a series of proton sponges, as well as to the related 1,8diaminonaphtalene and one of its derivatives that avoids the repulsion of the lone pairs of electrons. The results are presented in Table IV. The Jaguar pKa predictions performed without conformational search are inaccurate, again stresssing the importance of the necessity to perform conformational search for more stable results. The results obtained with conformational search are similar for all the protocols and show MUE’s of around 0.5 pKa units. As these compounds are ex-



Page 14 of 24 14

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

pected to posses very few conformations, it is not surprising that even the regular conformational search predicts pKa very well, and that an increase in the number of conformations does not lead to further improvement. This good performance aside, what is a truly remarkable achievement of the DFT-based method is that it reproduces the unusual experimental pKa trends perfectly. There is an enormous jump in pKa (7.49 units) when going from compound S1 to a similar but highly strained derivative S2. Jaguar pKa predicts this transition with almost quantitative accuracy (7.00 units using the 10-conformation gas phase optimization protocol, and 6.39 units using the 10-conformation solution phase optimization protocol). In contrast, the empirical pKa prediction program Epik89–91 fails to describe the very large experimental difference of 7.49 pKa units, giving instead a difference amounting to only 1.08 pKa units. This modest increase of pKa is presumably due to the simulated electrondonating effect of the methyl groups. If we assume that the electronic effect is described by Epik correctly, the rest of the experimental pKa difference is caused by the steric strain effects, which Epik, and likely all other empirical pKa predictors, cannot account for unless they are specifically trained, perhaps on a case by case basis, to do so. For example, the empirical pKa predictor Marvin92 shows results very similar to Epik’s, failing to recognize the influence of the steric repulsion. It is important to emphasize that Epik in itself is a highly accurate pKa predictor, as has been demonstrated in its application to diverse drug-like compounds39 . Therefore, Epik’s performance on proton sponges is atypical and only serves to highlight the advantages that DFT-based programs may have when predicting pKa of compounds displaying unusual sterical strain effects. An earlier version of Marvin, applied to typical drug-like structures, has also demonstrated average errors much lower than those shown in Table IV for proton sponges.93 Likewise, we conclude that Marvin’s performance on proton sponges does not reflect its typical accuracy on more common types of structures. Going back to the Table IV entries, there is also an extraordinary change in experimental pKa for two structurally very similar compounds, S2 and S3, which differ only by a CH2 group. The strain in S3, localized in the CH2 CH2 group, causes the lone pairs of electrons on the nitrogen atoms to avoid repulsion, so there is not as much energy gain upon protonation as in S2. Again, Jaguar pKa describes the experimental pKa drop correctly qualitatively and quantitatively, while Epik and Marvin are insensitive to this effect. The influence of the steric effects on pKa for the rest of the compounds from Table IV is similarly correctly reproduced by Jaguar pKa, but not by Epik and Marvin.

C.

Guanidines

The application of Jaguar pKa to a test set of guanidines from Ref. 94, apart from illustrating the use of Jaguar pKa computational protocols, elucidates a salient point about the use of the shell model in practice. Table V summarizes the results obtained for these guanidines.

Predictably, structures with minimal conformational flexibility (1a and 2j in Table V) are either not sensitive or only marginally sensitive to the number of conformations used in the calculations. More flexible structures always display some variation in the predicted pKa values, although all the structures from Table V contain very few rotatable bonds. Another observation is that the quality of the results depends strongly on the association of the functional group with a particular shell. Five shells match the R′′ N = C(NHR)(NHR′ ) functional group: “guanidine”, “amidine”, and “partly substituted amidine”, “protonation of sp2-like aliphatic nitrogen” and “protonation of generic atom”. The latter two shells are considered overly generic and are of little interest when more specific shells are matched. Therefore, we will only discuss the first three shells. The shell of “guanidine” produces the largest errors, the shell of “amidine” is more accurate, and the shell of “partly substituted amidine” is the most accurate, producing very low errors for multiconformational calculations that employ solution phase optimization. The problem is that the “guanidine” shell possesses the lowest RMSD for the linear fit of its parameters, and as such is currently automatically preferred by the Jaguar pKa program. The shell that produces the best results (“partly substituted amidine”) actually comes third in the ranking of the shells, because its parameter fitting has an RMSD (0.746) larger than the RMSD of the two other shells, “guanidine” (with RMSD = 0.320) and “amidine” (with RMSD = 0.662). Usually the shell that has the smallest linear fit RMSD is more specific, more relevant than the other shells to the structure under investigation, and produces the most accurate results. This is, however, not true for the test set from Table V. Upon reviewing our training set of guanidines used for generating the a, b parameters of the “guanidine” shell, we observed that the training set was not particularly representative of the structures from Table V, beyond both sets being guanidines. In particular, the training set, comprising 8 molecules, included only one structure with the ArN = C(NHR)(NHR′ ) moiety, with the rest of the structures having the sp2 nitrogen atom connected to an aliphatic carbon atom or hydrogen atom. From chemical intuition (compare pKa ’s of aliphatic and aromatic amines) we expect that the aromatic ring directly bonded to the the sp2 nitrogen atom would have a strong influence on the protonation properties of that atom. Our training sets for amidines and partly substituted amidines were even less representative of the guanidines from Table V. So the fact that the amidine shells were producing more accurate results was likely a coincidence. The problem comes down to the scarcity of experimental results that can be used for training the shell model. However, as the shell model is trained using more experimental data, with greater diversity of structures, we expect to encounter fewer cases where a less specific shell yields more accurate pKa ’s. It would be important to verify that there is nothing inherently wrong with the parameterization of our “guanidine” shell, and that the disappointing result produced with this shell in Table V was due to a bad representation in the training set. For this purpose, we constructed a small test of “ordinary” guanidines and applied Jaguar pKa to it. The results are as-


Page 15 of 24


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

sembled in Table VI. Each structure from Table VI matched at least two shells - “guanidine” and “amidine”. This should be sufficient for the intended verification, as the “guanidine” shell was performing significantly worse than the “amidine”, when applied to ArN = C(NHR)(NHR′ ) compounds. The results obtained with these guanidines, which should be closer in structure to those from the guanidine training set, do confirm that the “guanidine” shell performs clearly better than the “amidine” shell. These results corresponding to the best-performing shell generally agree with the previous conclusions. The molecules with low conformational flexibility (as G2 in this test set) have pKa ’s that are almost independent of the number of conformations considered in the computational protocol, although taking more conformations does gradually lower the error. The protocol that forgoes conformational search and performs geometry optimizations in vacuum gives a surprisingly good result, although the analogous protocol but with optimization in solution performs much worse. The multi-conformational protocols with the solution phase optimization give the best results.

VII. CONCLUSIONS AND OUTLOOK

This work discusses a computational method based on DFT calculations, aimed at accurate pKa predictions applicable to large, flexible organic molecules. We argue that while empirical methods that employ structure-activity relationships possess undeniable advantages of speed and robustness, methods that depend on quantum chemical energies should capture more physics of complex proton dissociation processes and be more suitable for describing equilibria involving steric and conformational effects. We propose a computational workflow implemented in the program Jaguar pKa which, in comparison with an earlier workflow, treats conformations better on several levels: an initial selection of conformations, a subsequent quantum mechanical optimization of their structures, and a statistical averaging of these conformations contributing to the observed pKa value. Our accounting for conformational entropy, and the theory behind its inclusion into micro-pKa ’s appears to be novel. As justification of the proposed treatment, we present applications for the new conformationally-aware approach to several sets of drug-like molecules. The results show clearly that a better treatment of conformations leads to more accurate results, with mean unsigned errors as low as 0.1-0.2 pKa units in some cases. Other pKa prediction protocols based on quantum chemical methods were able to show a comparable accuracy only for small molecules spanning an almost trivial conformational space14,15,17,37,95 or in application to a narrow chemical space.37 Errors obtained for more extended, conformationally flexible structures are typically much larger.7,23 Therefore we believe that the very high accuracy obtained in the present work for large, drug-like molecules with complex conformational landscapes is first of its kind and represents significant progress in the field of quantum chemical pKa prediction. Not only can very low errors be achieved, but we also demonstrate for the first time how the errors grow smaller with a

better conformational search conducted prior to quantum mechanical calculations and with a greater number of conformations used explicitly in quantum mechanical calculations. We recommend conducting conformational searches prior to performing a Jaguar pKa calculation, or using a builtin MacroModel conformational search protocol, even though this procedure requires additional computational resources. Failing to do so may have grave consequences for the accuracy of pKa prediction. As we have seen from the illustrative calculations, Jaguar pKa might happen to give very good pKa predictions even without a prior conformational search (Tables II, V, and VI). The results of such calculations are expected to depend heavily on the input conformations (see Figure 5) and are therefore not trustworthy. The use of multiple conformations consistently reduces either the the number of outliers, or both, with the exception of rigid or relatively rigid molecules, where explicit consideration of a single conformation is enough for good results. This is the main achievement of this work, showing how the quality of DFT-based pKa predictions can be improved systematically in challenging cases of drug-like flexible molecules. It cannot be firmly concluded from the present data which Jaguar pKa optimization protocols, those involving gas phase or solution phase optimizations, lead to more accurate results. For guanidines, solution phase optimizations performed clearly better (see Tables V, VI), but not for tertiary amines (see Table III), where gas phase optimizations gave more accurate results. The other test sets (Tables I, II, IV, and Figure 5) were either inconclusive or lacked data for proper comparison. Errors due to an inadequate treatment of conformational space are of course only one type of error that the quantum chemical pKa prediction may suffer from. Reducing this particular error, for example by such means as proposed in this work, is expected to improve the quality of pKa predictions across all chemical space, with the unimportant exception of small molecules where the errors are already low. However, other types of errors not related to the problem of treatment of conformations exist and are not uncommon to show up in practical applications. Errors due to insufficient specificity of the parameters a, b may present themselves when working with exotic structures or even with structures which do not share the essential structural elements with the structures from the relevant training set (as was illustrated by the results in Table V). Regular updating of the shell model with the newest literature data is therefore important. But even if such updates become frequent and the parameterization covers explicitly an ever growing fraction of the chemical space, there is still a problem of creating shells and assigning structures to the training sets, actions which are currently based on human, subjective decisions. Developing an algorithm that would create training sets, set up shells, and decide on the corresponding SMARTS patterns in a more objective manner will be a major improvement. Another type of error appears when acting on zwitterions. The existing shell model does contain shells and parameters specific to common zwitter-ionic structures, such as − NR+ 3 ...COO , but the coverage of the zwitter-ionic chemical space is not extensive, at least in part because of the short-



Page 16 of 24 16

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 I: The pKa values predicted by Jaguar pKa using different computational protocols. ‘No CS’ stands for no conformational search performed, with the calculation being conducted on the input conformation. In ‘CS’, regular conformational search was performed on the input conformation. In ‘N conf’, an aggressive CS was performed, N (e.g., N = 1, 10) conformations were selected, and a geometry optimization in solution phase was conducted. The structure indices follow the numbering from Ref. 78. The predicted pKa is considered an outlier if its difference in value is equal to or larger than 1.0 pKa unit, in comparison with the experimental measurement. 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 R2 − R5 CN CN R3 = F CN R2 = F CN CN R4 = CF3 CN R5 = CF3 CN R3 = F CN R5 = CF3 CN R5 = CF3 , R3 = F CN R4 = CF3 , R3 = F CN R2 = F, R3 = F CN R2 = F, R3 = F Cl Cl R2 = F, R3 = F Cl Cl R2 = F, R3 = F MUE Outliers

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

R5

X Exp. pKa No CS O 9.7 9.1 O 8.1 6.6 O 7.4 6.0 O 7.3 7.3 O 7.3 6.6 O 7.0 6.8 O 6.7 4.4 O 6.3 4.7 O 5.9 5.5 O 5.8 3.8 O 5.8 5.1 O 5.1 3.2 O 9.8 9.2 O 6.3 5.1 S 9.0 7.8 S 6.1 5.3 1.07 8

age of the experimental data that could be used for parameterization. Our tests show that the parameters from the generic shell applied to zwitter-ions not explicitly covered by the shell model yield errors substantially larger than those observed for non-zwitter-ionic structures. Both the parameters of the implicit solvation model in application to zwitter-ions and the corresponding parameters of the shell model in the pKa protocol might need to be revisited and improved. Errors arising due to multiple possible tautomers for the given structure are culprits behind a number of outliers: either the user fails to act on the dominant conformer, or more than one conformer co-exist in solution and contribute to observable pKa with the user operating Jaguar pKa maybe aware of only one conformer. Automatically generating, scoring, and acting on different tautomeric forms in one pKa prediction workflow will be another major improvement. Finally, there is a similar problem of multiple protonation channels. In a molecule containing several protonatable or deprotonatable functional groups, it is not always clear a priori, even for the experts, which of those groups is or are responsible for the observable pKa value. Significant mispredictions may result when the user applies Jaguar pKa to an irrelevant protonation channel or fails to account for all relevant protonation channels.

CS 1 conf. 10 conf. 9.0 9.4 9.7 7.7 8.1 8.2 6.0 7.6 7.4 7.3 7.4 7.7 5.4 6.3 7.0 6.8 7.4 7.3 5.0 7.2 7.3 6.1 6.2 6.3 5.8 5.9 6.0 5.2 5.9 6.0 5.8 6.5 6.3 4.3 5.5 4.7 9.4 9.5 9.8 6.2 6.6 6.4 9.5 8.9 9.0 6.1 6.2 6.3 0.56 0.29 0.20 3 1 0

The number and scope of the remaining problems associated with quantum chemical pKa prediction suggest that much more scientific and engineering work is needed before very accurate and robust pKa prediction covering all practically important aspects of protonation processes becomes available. We would like to believe that the present work makes an important step in the direction of that necessary progress.

VIII. ACKNOWLEDGMENT

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

IX. SUPPLEMENTARY INFORMATION

The hierarchical organization of the current shell model used for classification of the functional groups as well as average errors of parameterization produced by different parameterization models are available in the supplementary information.


Page 17 of 24


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 II: Standard deviations computed for the set of 15 pKa values obtained in 15 separate Jaguar pKa calculations that used different and distinct input conformations. The RMSD range was 0.92 to 4.37 Angstrom for all relevant pairs. The ideal result would be zero deviation. Results for three molecules from Ref. 80 using different Jaguar pKa protocols are presented. The total number of outliers with each computational protocol is determined for the set of 45 pKa predictions. An outlier is defined as a pKa prediction deviating from the experimental value by 1.0 pKa units or more.

Structure (2S,3R)-7a (2R,3R)-7a 3R-6a Outliers

Gas phase optimization No CS Reg. CS 1 conf. 5 conf. 2.092 2.011 0.423 0.091 0.692 0.853 0.052 0.048 1.770 0.343 0.024 0.032 9 8 0 0

10 conf. 0.099 0.041 0.034 0

Solution phase optimization No CS 1 conf. 5 conf. 10 conf. 0.972 0.171 0.049 0.056 0.609 0.563 0.054 0.067 0.519 0.159 0.051 0.047 8 1 0 0



Page 18 of 24 18

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 III: pKa values predicted with various Jaguar pKa computational protocols on a test set of tertiary amines extracted from Ref. 82. Mean unsigned errors (MUE) and the number of outliers are also reported for each protocol. An outlier is defined as a pKa prediction deviating from the experimental value by 1.0 pKa units or more. O

N

A1

N

O

A8 B1

O N N

A2

A9 B2

HO O N

N

A3

A10

B3 O

HO O

N

A4

N

A11

N

A12

n

H2 C

B4

HO N

N

A5

O

N

F

O

N

A6 O

O

N

A7

N

A13

N

A14

O

n B5

N N

N

B6 N

R

B7

O

Index 76 78 46 51 79 80 81 82 83 84 85 86 87 88 89 90 91 94 95 96 97 98 99 100 101 103 104

Structure A B n A9 B2 0 A8 B7 0 A9 B7 1 A9 B6 1 A8 B5 0 A8 B4 0 A8 B3 0 A5 B2 1 A5 B5 0 A5 B5 0 A6 B2 1 A6 B1 0 A11 B2 1 A11 B2 1 A10 B2 0 A10 B2 0 A12 B2 1 A13 B2 0 A2 B2 0 A2 B2 0 A7 B2 0 A14 B2 0 A1 B2 0 A1 B2 0 A3 B2 0 A4 B2 0 A4 B2 0 MUE Outliers

Gas phase optimization R No CS Reg. CS 1 conf. 5 conf. 10 conf. OCH3 9.43 8.13 8.98 9.00 9.11 OCH3 7.94 8.23 7.99 7.85 7.90 OCH3 8.23 8.68 8.69 8.58 8.69 OCH3 9.02 9.30 9.17 9.01 9.13 OCH3 7.81 9.04 8.44 8.23 8.31 OCH3 8.49 8.34 8.44 8.49 9.43 H 8.68 8.62 8.64 8.72 8.66 OCH3 6.60 7.96 7.99 8.29 8.18 OCH3 7.34 7.03 7.49 7.50 7.79 H 7.80 7.22 7.59 7.86 7.78 OCH3 7.96 7.91 7.75 8.39 8.25 OCH3 8.50 8.04 7.81 7.88 7.89 H 8.33 8.12 8.60 8.48 8.61 OCH3 8.76 9.05 8.82 8.29 8.83 H 8.72 8.33 8.85 8.65 8.62 OCH3 8.61 9.17 8.38 8.32 8.32 OCH3 8.24 9.41 9.45 9.30 9.27 OCH3 8.20 8.42 8.53 7.95 8.00 OCH3 8.15 7.69 3.87 7.00 7.11 H 6.54 8.22 6.66 7.36 7.71 OCH3 8.56 8.21 7.56 7.75 6.47 OCH3 7.78 7.07 6.98 6.94 6.93 OCH3 8.51 8.32 9.10 8.49 8.53 H 8.36 8.57 8.64 8.39 8.38 OCH3 7.98 7.88 9.58 8.55 8.35 OCH3 7.83 8.00 8.13 7.93 7.96 H 7.67 7.80 7.86 7.89 7.90 0.69 0.56 0.69 0.56 0.51 6 5 4 3 2

Solution phase optimization No CS 1 conf. 5 conf. 10 conf. Exp. pKa 9.38 9.08 8.73 8.70 9.9 6.19 7.40 7.48 7.50 8.5 8.72 8.69 8.15 8.47 9.2 9.65 9.13 9.43 9.33 9.1 8.10 8.44 7.81 7.89 8.6 9.17 8.44 8.51 8.49 9.3 5.20 8.64 8.49 8.64 9.3 8.34 7.99 7.62 7.68 8.5 6.85 7.49 7.06 7.41 8.4 6.95 7.11 7.42 7.48 8.4 7.60 7.37 7.58 7.58 6.5 8.05 8.63 8.21 8.20 7.5 9.42 8.60 8.21 7.67 8.9 8.49 8.82 7.67 7.95 8.9 8.77 8.85 7.95 8.68 9.0 8.17 7.87 8.68 8.05 9.0 8.11 9.45 8.60 9.60 9.0 7.50 8.53 7.93 7.93 8.5 8.44 3.87 7.83 7.83 8.0 7.68 7.02 7.70 7.70 8.5 8.21 7.56 7.58 7.58 8.1 6.92 5.18 6.17 6.17 6.0 7.68 9.10 8.05 8.05 8.6 6.69 8.27 7.84 7.84 8.6 7.74 9.58 7.76 7.36 8.4 7.29 8.56 7.79 7.75 8.2 6.98 7.86 7.75 7.73 8.0 0.93 0.76 0.73 0.70 8 7 6 6


Page 19 of 24


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 IV: pKa values predicted with various Jaguar pKa computational protocols on a test set of proton sponges. Mean unsigned errors (MUE) and the number of outliers are also reported for each protocol. An outlier is defined as a pKa prediction deviating from the experimental value by 1.0 pKa units or more. All experimental pKa data refer to measurements in water at room temperature. NH2

NH2

N

S1

N

N

S2

S3

O

N

N

S4

Index S1 S2 S3 S4 S5 S6 MUE Outliers

Gas phase optimization No CS Reg. CS 1 conf. 5 conf. 10 conf. 2.01 4.11 4.28 4.28 4.28 11.19 11.28 11.31 11.28 11.28 7.41 4.42 4.18 4.06 4.06 8.95 9.83 9.80 9.83 9.83 8.44 8.43 8.46 8.54 8.54 10.03 10.05 10.08 10.07 10.07 1.44 0.50 0.52 0.55 0.55 3 0 0 1 1

N

O

N

N

N

S5

N

S6

Solution phase optimization No CS 1 conf. 5 conf. 10 conf. 1.90 4.86 4.86 4.86 11.22 11.19 11.25 11.25 7.64 4.39 4.32 4.32 8.82 9.43 9.45 9.46 8.36 8.41 8.50 8.49 4.28 10.02 9.88 9.88 2.44 0.53 0.55 0.55 4 0 1 1

a Ref.

96 97 98 d Ref. 99 b Ref. c Ref.


Epik Marvin Exp. pKa 4.46 4.25 4.61a 5.54 5.01 12.10b 5.84 4.37 4.62 c 6.02 5.44 10.27c 1.86 1.49 7.49d 4.60 5.16 9.97c 3.86 3.61 5 4


Page 20 of 24 20

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 V: pKa values predicted with various Jaguar pKa computational protocols on a test set of guanidine structures extracted from Ref. 82. Mean unsigned errors (MUE) and the number of outliers are also reported for each protocol. An outlier is defined as a pKa prediction deviating from the experimental value by 1.0 pKa units or more. The predicted pKa values and their statistical summaries in plain font correspond to the “guanidine” shell of the shell model, while the figures in italics are associated with the “amidine” shell of the shell model, and the figures in bold font refer to the “partly substituted amidine” shell of the shell model. All experimental pKa data refer to measurements in water at room temperature. R1

CH3

R2 NH R4 N

N H

R3

Structure R2 R3

Index

R1

1a

Cl

H

H

H

1b

Cl

H

H

CH2 CHF2

2j

H

H

H

CH3

2o

H

H

H

CH2 CHF2

2p

H

H

H

CH2 CF3

3

Cl

Cl

H

CH2 CHF2

4

H

H OCH3 CH2 CHF2

6

7

OCH3 Cl

Cl

H

R4

H

CH2 CHF2

Cl

CH2 CHF2

MUE

Outliers

Gas phase optimization No CS Reg. CS 1 conf. 5 conf. 8.62 8.70 8.70 8.70 9.36 9.43 9.42 9.42 9.59 9.65 9.65 9.65 7.02 6.93 6.95 6.85 8.13 8.06 8.08 8.03 8.52 8.46 8.48 8.45 10.47 10.58 10.75 10.49 10.78 10.87 11.00 10.80 10.84 10.91 11.02 10.85 8.42 7.98 7.98 8.04 9.21 8.87 8.87 8.88 9.46 9.17 9.17 9.16 7.36 6.65 6.69 7.64 8.39 7.85 7.88 8.58 8.75 8.27 8.30 8.89 6.17 6.23 6.21 6.09 7.48 7.53 7.51 7.45 7.95 7.99 7.98 7.94 9.53 9.07 9.10 8.99 10.06 9.71 9.73 9.67 10.20 9.90 9.92 9.88 6.68 7.19 6.56 7.41 7.87 8.27 7.78 8.38 8.30 8.64 8.22 8.70 5.00 5.26 5.29 6.00 6.58 6.78 6.81 7.33 7.17 7.35 7.37 7.81 1.60 1.68 1.75 1.50 0.69 0.76 0.82 0.62 0.38 0.44 0.49 0.32 7 8 8 8 3 2 2 1 0 0 0 0

10 conf. 8.70 9.42 9.65 7.43 8.45 8.80 10.49 10.80 10.85 8.61 9.34 9.56 7.65 8.57 8.88 6.80 7.93 8.33 9.23 9.87 10.06 6.65 7.84 8.27 6.11 7.41 7.89 1.34 0.50 0.23 7 1 0

Solution phase optimization No CS 1 conf. 5 conf. 10 conf. Exp. pKa 8.77 8.80 8.80 8.80 9.48 9.50 9.50 9.50 9.9 9.70 9.72 9.72 9.72 6.01 8.01 7.54 7.69 7.36 8.89 8.52 8.68 8.9 7.85 9.19 8.86 9.01 10.50 10.61 10.28 10.28 10.81 10.89 10.63 10.63 10.6 10.86 10.93 10.70 10.70 7.73 9.37 8.75 8.86 8.68 9.94 9.47 9.58 9.7 9.00 10.10 9.69 9.80 7.75 7.97 7.93 7.84 8.69 8.86 8.84 8.76 9.2 9.01 9.16 9.14 9.08 5.00 7.35 6.93 7.09 6.58 8.39 8.04 8.20 8.5 7.17 8.74 8.44 8.59 8.23 9.11 9.09 9.15 9.06 9.74 9.75 9.82 10.2 9.34 9.92 9.95 10.03 7.46 7.81 7.34 7.68 8.47 8.47 8.40 8.60 8.9 8.82 9.05 8.76 8.91 3.26 5.72 6.10 6.13 5.25 7.14 7.39 7.45 7.8 6.01 7.65 7.86 7.93 2.11 1.00 1.22 1.13 1.08 0.30 0.36 0.28 0.72 0.23 0.10 0.11 8 6 7 7 5 0 0 0 3 0 0 0


Page 21 of 24


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 VI: pKa values predicted with various Jaguar pKa computational protocols on a test set of guanidine structures of the type X − N = C(NHR)(NHR′ ), where X is non-aryl. Mean unsigned errors (MUE) and the number of outliers are also reported for each protocol. An outlier is defined as a pKa prediction deviating from the experimental value by 1.0 pKa units or more. The predicted pKa values and their statistical summaries in plain font correspond to the “guanidine” shell of the shell model, while the figures in italics are associated with the “amidine” shell of the shell model. All experimental pKa data refer to measurements in water at room temperature. O

N

N

NH2 N NH

N

O

NH2

G1 OH H

G2 H

HO

OH NH2

O

N

OH H

H2N

H N N H

O OH

N HO

H

O

H

G3 Gas phase optimization Index No CS Reg. CS 1 conf. 5 conf. G1 11.39 11.15 11.06 11.15 11.42 11.23 11.17 11.24 G2 4.93 6.08 6.16 6.16 6.53 7.41 7.47 7.47 G3 9.11 9.86 9.76 10.59 9.74 10.31 10.24 10.88 G4 9.75 10.76 10.96 9.45 10.09 10.86 11.02 9.88 MUE 0.66 0.92 0.94 0.79 0.90 1.41 1.43 1.32 Outliers 1 2 2 1 1 3 3 3

NH

G4

10 conf. 11.10 11.21 6.16 7.47 10.61 10.89 9.48 9.92 0.79 1.33 1 3

Solution phase optimization No CS 1 conf. 5 conf. 10 conf. Exp. pKa 9.88 10.34 10.34 10.33 11.00a 10.26 10.61 10.61 10.58 6.77 5.93 5.93 5.93 5.75b 7.94 7.30 7.30 7.30 10.43 9.22 9.25 9.23 8.76c 10.75 9.82 9.87 9.86 10.62 10.60 9.81 9.78 8.68d 12.27 10.74 10.18 10.16 1.44 0.81 0.62 0.61 2.13 1.27 1.14 1.14 4 1 1 1 3 3 3 3

a Ref.

100 101 102 d Ref. 103 b Ref. c Ref.



Page 22 of 24 22

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

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

Electronic address: art.bochevarov@schrodinger. com Alongi, K. S.; Shields, G. C. Ann. Rep. Comp. Chem. 2010, 6, 113–138. Ho, J.; Coote, M. L. Theor. Chem. Acc. 2010, 125, 3–21. Eckert, F.; Diedenhofen, M.; Klamt, A. Mol. Phys. 2010, 3-4, 229–241. Liao, C.; Niklaus, M. C. J. Chem. Inf. Model. 2009, 49, 2801– 2812. Ho, J. Aus. J. Chem. 2014, 67, 1441–1460. Settimo, L.; Bellman, K.; Knegtel, R. M. A. Pharm. Res. 2014, 114, 1082–1095. Kaupmees, K.; Trummal, A.; Leitoa, I. Croat. Chem. Acta 2014, 87, 385–395. Rajapakse, H. A. et al. Bioorg. Med. Chem. Lett. 2010, 20, 1885– 1889. Lanevskij, K.; Japertas, P.; Didziapetris, R.; Petrauskas, A. J. Pharm. Sci. 2009, 98, 122–134. Sprous, D. G.; Palmer, R. K.; Swanson, J. T.; Lawless, M. Curr. Top. Med. Chem. 2010, 10, 619–637. Cheng, J.; Sprik, M. J. Chem. Theory Comput. 2010, 6, 880–889. Bryantsev, V. S. Chem. Phys. Lett. 2013, 558, 42–47. Gallus, D. R.; Wagner, R.; Wiemers-Meyer, S.; Winter, M.; Cekic-Laskovic, I. Electrochim. Acta 2015, 184, 410–416. Liptak, M. D.; Shields, G. C. J. Am. Chem. Soc. 2001, 123, 7314– 7319. Liptak, M. D.; Gross, K. C.; Seybold, P. G.; Feldgus, S.; Shields, G. C. J. Am. Chem. Soc. 2002, 124, 6421–6427. Klicić, J. J.; Friesner, R. A.; Liu, S.-Y.; Guida, W. C. J. Phys. Chem. A 2002, 106, 1327–1335. Casasnovas, R.; Fernández, D.; Ortega-Castro, J.; Frau, J.; Donoso, J.; Muñoz, F. Theor. Chem. Acc. 2011, 130, 1–13. Jerome, S. V.; Hughes, T. F.; Friesner, R. F. J. Phys. Chem. B 2014, 118, 80088016. Lee, A. C.; yu Yu, J.; Crippen, G. M. J. Chem. Inf. Model. 2008, 48, 2042–2053. Fraczkiewicz, R.; Lobell, M.; Göller, A. H.; Krenz, U.; Schoenneis, R.; Clark, R. D.; Hillisch, A. J. Chem. Inf. Model. 2015, 55, 389–397. Rebollar-Zepeda, A. M.; Galano, A. Int. J. Quant. Chem 2012, 112, 3449–3460. Saracino, G. A. A.; Improta, R.; Barone, V. Chem. Phys. Lett. 2003, 373, 411–415. Ho, J.; Coote, M. L.; Franco-Pérez, M.; Gómez-Balderas, R. J. Phys. Chem. A 2010, 44, 11992–12003. Harding, A. P.; Popelier, P. L. A. Chem. Phys. Phys. Chem. 2011, 23, 11264–11282. Harding, A. P.; Popelier, P. L. A. Chem. Phys. Phys. Chem. 2011, 23, 11283–11293. Griffiths, M. Z.; Alkorta, I.; Popelier, P. L. A. Mol. Inf. 2013, 32, 363–376. Ugur, I.; Marion, A.; Parant, S.; Jensen, J. H.; Monard, G. J. Chem. Inf. Model. 2014, 54, 22002213. Yu, H.; Wondrousch, D.; Yuan, Q.; Lin, H.; Chen, J.; Hong, H.; Schüu¨ rmann, G. Chemosphere 2015, 138, 829–836. Sastre, S.; Casasnovas, R.; Muñoz, F.; Frau, J. Phys. Chem. Chem. Phys. 2016, 18, 11202–11212. Crugeiras, J.; R´ıos, A.; Maskill, H. J. Phys. Chem. A 2011, 115, 12357–12363. Feng, S.; Bagia, C.; Mpourmpakis, G. J. Phys. Chem. A 2013, 117, 5211–5219.

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

Alexov, E.; Mehler, E. L.; Baker, N.; Baptista, A. M.; Huang, Y.; Milletti, F.; Nielsen, J. E.; Farrell, D.; Carstensen, T.; Olsson, M. H. M.; Shen, J. K.; Warwicker, J.; Williams, S.; Word, J. M. Proteins 2011, 79, 3260–3275. Andersson, M. P.; Jensen, J. H.; Stipp, S. L. S. PeerJ 2013, 1, e198. Stanton, C. L.; Houk, K. N. J. Chem. Theory Comput. 2008, 4, 951–966. Meyer, T.; Knapp, E.-W. J. Chem. Theory Comput. 2015, 11, 2827–2840. Koch, W.; Holthausen, M. C. A Chemist’s Guide to Density Functional Theory; Wiley-VCH Verlag: Weinheim, Germany, 2002. Bryantsev, V. S.; Diallo, M. S.; Goddard III, W. A. J. Phys. Chem. A 2007, 111, 4422–4430. Morgenthaler, M.; Schweizer, E.; Hoffmann-Röder, A.; Benini, F.; Martin, R. E.; Jaeschke, G.; Wagner, B.; Fischer, H.; Bendels, S.; Zimmerli, D.; Schneider, J.; Diederich, F.; Kansy, M.; Müller, K. Chem. Med. Chem. 2007, 6, 1100–1115. Settimo, L.; Bellman, K.; Knegtel, R. M. Pharm. Res. 2014, 31, 1082–1095. 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. Int. J. Quant. Chem. 2013, 113, 2110– 2142. Schrödinger Release 2016-2: Jaguar, version 9.2. Schrödinger, LLC: New York, NY, 2016. Lim, C.; Bashford, D.; Karplus, M. J. Phys. Chem. 1991, 95, 5610–5620. Jang, Y. H.; Sowers, L. C.; C ¸ aˇgin, T.; Goddard III, W. A. J. Phys. Chem. A 2001, 105, 274–280. Tissandier, M. D.; Cowen, K. A.; Feng, W. Y.; Gundlach, E.; Cohen, M. H.; Earhart, A. D.; Coe, J. V.; Tuttle, Jr., T. R. J. Phys. Chem. A 1998, 102, 7787–7794. Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2006, 110, 16066–16081. Cortis, C. M.; Friesner, R. A. J. Comp. Chem. 1997, 18, 1570– 1590. Cortis, C. M.; Friesner, R. A. J. Comp. Chem. 1997, 18, 1591– 1608. Marten, B.; Kim, K.; Cortis, C.; Friesner, R. A.; Murphy, R. B.; Ringnalda, M. N.; Sitkoff, D.; Honig, B. J. Phys. Chem. 1996, 100, 11775–11788. Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215–241. Silva, P. J.; Ramos, M. J. Comp. Theor. Chem. 2011, 966, 120– 126. ˇ acˇ , J.; Hobza, P. Valdes, H.; Pluhácˇ ková, K.; Pitonák, M.; Rez´ Phys. Chem. Chem. Phys. 2008, 10, 2747–2757. Zhang, I. Y.; Wu, J.; Luo, Y.; Xu, X. J. Comp. Chem. 2011, 9, 1824–1838. Johnson, E. R.; Clarkin, O. J.; DiLabio, G. A. J. Phys. Chem. A 2003, 107, 9953–9963. Chandra, A. K.; Uchimaru, T. Int. J. Mol. Sci. 2002, 3, 407–422. Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. J. Chem. Theory Comput. 2013, 9, 609–620. Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. A 2006, 110, 2493–2499. Marenich, A. V.; Ding, W.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. Lett. 2012, 3, 1437–1442. Casasnovas, R.; Ortega-Castro, J.; Donoso, J.; Frau, J.; Muñoz, F. Phys. Chem. Chem. Phys. 2013, 38, 16303–16313. Kromann, J. C.; Larsen, F.; Moustafa, H.; Jensen, J. H. PeerJ


Page 23 of 24


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

60

61

62

63

64

65

66

67

68

69

70 71

72

73

74

75

76 77

78 79

80 81 82 83

2016, 8, e2335. Churakov, S. V.; Labbez, C.; Pegado, L.; Sulpizi, M. J. Phys. Chem. C 2014, 118, 11752–11762. Uddin, N.; Choi, T. H.; Choi, C. H. J. Phys. Chem. B 2013, 117, 6269–6275. Miguel, E. L. M.; Silva, P. L.; Pliego, J. R. J. Phys. Chem. B 2014, 118, 5730–5739. Jackson, V. E.; Felmy, A. R.; Dixon, D. A. J. Phys. Chem. A 2015, 119, 2926–2939. Mari´ın-Luna, M.; Alkorta, I.; Elguero, J. New J. Chem. 2015, 39, 2861–2871. Xue, X.-S.; Wang, Y.; Yang, C.; Ji, P.; Cheng, J.-P. J. Org. Chem. 2015, 80, 8997–9006. Klamt, A.; Eckert, F.; Diedenhofen, M.; Beck, M. E. J. Phys. Chem. A 2003, 107, 9380–9386. Brown, T. N.; Mora-Diez, N. J. Phys. Chem. B 2006, 110, 9270– 9279. Farràs, P.; Teixidor, F.; Branchadell, V. Inorg. Chem. 2006, 45, 7947–7954. Casasnovas, R.; Frau, J.; Ortega-Castro, J.; Salvà, A.; Donoso, J.; Muñoz, F. Int. J. Quant. Chem. 2010, 110, 323–330. Rossini, E.; Knapp, E.-W. J. Comp. Chem. 2016, 12, 1082–1091. James, C. A.; Weininger, D. Daylight Theory Manual, Version 4.9. Daylight Chemical Information Systems, Inc.: Laguna Niguel, CA, 2011. Schrödinger Release 2016-1: MacroModel, version 11.2. Schrödinger, LLC: New York, NY, 2016. Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657–1666. Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225–11236. Shivakumar, D.; Williams, J.; Wu, Y.; Damm, W.; Shelley, J.; Sherman, W. J. Chem. Theory Comput. 2010, 6, 1509–1519. Harder, E. et al. J. Chem. Theor. Comput. 2016, 12, 281–296. Csizmadia, F.; Tsantili-Kakoulidou, A.; Panderi, I.; Darvas, F. J. Pharm. Sci. 1997, 86, 865–871. Hilpert, H. et al. J. Med. Chem. 2013, 56, 3980–3995. Tresadern, G.; Delgado, F.; Delgado, O.; Gijsen, H.; Macdonald, G. J.; Moechars, D.; Rombouts, F.; Alexander, R.; Spurlino, J.; Gool, M. V.; Vega, J. A.; Trabanco, A. A. Bioorg. Med. Chem. Lett. 2011, 21, 7255–7260. Rombouts, F. J. R. et al. J. Med. Chem. 2015, 58, 8216–8235. Ginman, T. et al. J. Med. Chem. 2013, 56, 4181–4205. Johansson, A. et al. J. Med. Chem. 2016, 59, 2497–2511. Alder, R. W.; Bowman, P. S.; Steele, W. R. S.; Winterman, D. R.

84 85

86

87

88

89

90

91

92

93

94

95

96 97

98

99

100

101

102

103

Chem. Comm. 1968, 723–724. Alder, R. W. Chem. Rev. 1989, 89, 1215–1223. Llamas-Saiz, A. L.; Foces-Foces, C.; Elguero, J. J. Mol. Struct. 1994, 328, 297–323. Margetić, D.; Ishikawa, T.; Kumamoto, T. Eur. J. Org. Chem. 2010, 2010, 6563–6572. Ozeryanskii, V. A.; Gorbacheva, A. Y.; Pozharskii, A. F.; Vlasenko, M. P.; Tereznikov, A. Y.; Chernov’yants, M. S. Org. Biomol. Chem. 2015, 13, 8524–8532. Belding, L.; Stoyanov, P.; Dudding, T. J. Org. Chem. 2010, 81, 6–13. Schrödinger Release 2016-2: Epik, version 3.6. Schrödinger, LLC: New York, NY, 2016. Shelley, J. C.; Cholleti, A.; Frye, L.; Greenwood, J. R.; Timlin, M. R.; Uchimaya, M. J. Comput. Aided Mol. Des. 2010, 21, 681–691. Greenwood, J. R.; Calkins, D.; Sullivan, A. P.; Shelley, J. C. J. Comput. Aided Mol. Des. 2010, 24, 591–604. Marvin: 16.7.18.0. ChemAxon (http://www.chemaxon.com), 2016. Balogh, G. T.; Tarcsay, A.; Keser˝u, G. M. J. Pharm. Bio. Anal. 2012, 67-68, 63–70. Peters, J.-U.; Lübbers, T.; Alanine, A.; Kolczewski, S.; Blasco, F.; Steward, L. Bioorg. Med. Chem. Lett. 2008, 18, 262– 266. Trummal, A.; Rummel, A.; Lippmaa, E.; Koppel, I.; Koppel, I. A. J. Phys. Chem. A 2011, 115, 6641–6645. Alexander, A. J.; Kebarle, P. Anal. Chem. 1986, 58, 471–478. Hibbert, F.; Kenneth, P. P. J. Chem. Soc., Perkin Trans. 2 1983, 1895–1900. Hibbert, F.; Simpson, G. R. J. Chem. Soc., Perkin Trans. 2 1987, 613–616. Hibbert, F.; Simpson, G. R. J. Chem. Soc., Perkin Trans. 2 1987, 243–246. Maillard, M. C.; Perlman, M. E.; Amitay, O.; Baxter, D.; Berlove, D.; Connaughton, S.; Fischer, J. B.; Guo, J. Q.; Hu, L.Y.; McBurney, R. N.; Nagy, P. I.; Subbarao, K.; Yost, E. A.; Zhang, L.; Durant, G. J. J. Med. Chem. 2011, 41, 3048–3061. Chernyshev, V. M.; Chernysheva, A. V.; Starikova, Z. A. Heterocycles 2010, 81, 2291–2311. Goto, T.; Kishi, Y.; Takahashi, S.; Hirata, Y. Tetrahedron 1965, 21, 2059–2088. Tarasova, E. V.; Chernyshev, V. M.; Chernysheva, A. V.; Abagyan, R. S. Russ. J. Appl. Chem. 2011, 84, 400–406.



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

137x148mm (144 x 144 DPI)


Page 24 of 24

Multiconformation, Density Functional Theory ... - ACS Publications

Recommend Documents