Reliable and Performant Identification of Low-energy Conformers in

Subscriber access provided by University of Winnipeg Library

Computational Chemistry

Reliable and Performant Identification of Lowenergy Conformers in Gas-phase and Water Anna Theresa Cavasin, Alexander Hillisch, Felix Uellendahl, Sebastian Schneckener, and Andreas H. Göller J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00151 • Publication Date (Web): 02 May 2018 Downloaded from http://pubs.acs.org on May 10, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a 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 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.

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

Journal of Chemical Information and Modeling

Reliable and Performant Identification of Lowenergy Conformers in Gas-phase and Water Anna Theresa Cavasin1 $, Alexander Hillisch1, Felix Uellendahl1 #, Sebastian Schneckener2, Andreas H. Göller1* 1

Bayer AG, Drug Discovery, Chemical Research, 42096 Wuppertal, Germany

2

Bayer AG, Engineering & Technology, Applied Mathematics, 51368 Leverkusen, Germany

ABSTRACT Prediction of compound properties from structure via QSAR and machine-learning approaches is an important computational chemistry task in small molecule drug research. Though many such properties are dependent on 3D structures or even conformer ensembles, the majority of models are based on descriptors derived from 2D structure. Here, we present results from a thorough benchmark study of force field, semiempirical and density functional methods for the calculation of conformer energies in gas-phase and water solvation as a foundation for the correct identification of relevant low-energy conformers. We find that the tight-binding ansatz GFN-xTB shows the lowest error metrics and highest correlation to the benchmark PBE0D3(BJ)/def2-TZVP in gas-phase for the computationally fast methods, and that OPLS3 in solvent becomes comparable in performance. MMFF94, AM1, and DFTB+ perform worse, whereas the performance-optimized but far more expensive functional PBEh-3c yields energies

ACS Paragon Plus Environment

1

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

Page 2 of 49

almost perfectly correlated to the benchmark and should be used whenever affordable. Based on our findings, we have implemented a reliable and fast protocol for the identification of lowenergy conformers of drug like molecules in water which can be used for the quantification of strain energy and entropy contributions to target binding as well as for derivation of conformerensemble dependent molecular descriptors.

INTRODUCTION Many properties and experimental measures of the chemical structures relevant in pharmaceutical and agrochemical research are dependent on low-energy conformer ensembles. Whereas structure-based design and pharmacophore modeling is done in 3D space, multiple properties in the area of ADMET (absorption, distribution, metabolism, excretion, toxicology), e.g. solubility or permeation, are predicted by 2D QSAR, therefore ignoring conformational aspects completely. We and others1,2 expect that this simplification has its limitations and that for next generation ADMET prediction algorithms 3D conformers or conformer ensembles will play some role, especially for the so-called beyond rule of five (bRo5) compounds.3 Some first steps in the direction of 3D ADMET are approaches like physics-based permeation modeling.4,5 The main reason for the limited utilization of 3D structures is the lack of a quick and reliable process for the energy calculation of generated conformer sets. With the advent of quantum chemistry-derived accurate force fields for small molecules, the 3D conformer problem can be readdressed. We have recently published a first step into 3D ADMET. Via machine-learning with chargedistribution descriptors from gas-phase 3D structures6 obtained from one single conformer we


2

Page 3 of 49 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


have implemented a robust and performant method for the site-of-metabolism prediction of cytochrome P450 catalyzed transformations.7 A general understanding of the regioselectivity of chemical reactions will nevertheless often require the low-energy conformer ensemble in solvent. We are currently exploring possibilities to apply descriptors derived from ensembles for the prediction of ADMET properties where flexibility and the surrounding medium plays a major role, such as membrane permeation or solubility. Indeed, the surrounding medium is a major hurdle in this respect. Conformer structures can be determined experimentally in any aggregate phase and medium but experiments in different phases or media will give different results, for well-understood reasons. Gas-phase coordinates describe the structure in a more or less undisturbed state; the experiments are limited to relatively small structures which can sublimate into gas-phase without decomposition; solvent-phase coordinates can only indirectly be determined by shifts and couplings from spectroscopic methods like NMR or IR. Solid state coordinates are obtained by crystallography either for the ligand itself or for a ligand cocrystallized with a target protein. Small molecule crystals provide high-resolution coordinates which however often do not represent the global minimum conformation, as they are defined by intra - and more importantly - intermolecular interactions like hydrogen bonds, pi-stacking, dispersion, charge-charge interactions etc. which strongly influence the torsional angles in particular. Coordinates derived from protein-ligand complexes on the other hand are significantly less precise and accurate, providing only heavy-atom positions which often have non-equilibrium distances, angles and torsions,8 and even high-resolution structures often have no electron density for parts of the ligand.9 Additionally, a study by Perola10 reported that from the 150 protein-ligand complexes evaluated, about 60% were no local minima, about 60% had strain energies of up to 5 kcal/mol and at least 10% had strain energies higher than 9 kcal/mol. Other


3


Page 4 of 49

studies, using higher levels of theory, report much lower (10 kcal/mol) strain energies, as summarized by Hawkins.1111 Regardless of whether researchers work in the pharmaceutical industry or academia they should aim to design compounds which are able to optimally adopt the bio-active conformation readily achieved according to Perola’s work - and at the same time are not too flexible. Both aspects avoid loss of binding free energy due to strain energy12 and entropy.13 On the other hand, molecules should not be too flat14 since they have to be soluble. The problem we are dealing with, from a computational chemist’s standpoint, consists of two steps. First, raw conformer ensembles have to be created that ideally capture the complete accessible low-to-medium energy conformer space (energy thresholds between 3 and 10 kcal/mol are commonly applied15,16,17,18). We aim at a low overall number of representative, yet diverse conformers since all have to be geometry-optimized in the second step to derive the final coordinates and energies. Since we commonly deal with sets upto thousands of molecules and work under severe time-constraints, the overall process has to be robust, and both compute- and wall-time efficient. Concerning step one, there are multiple publications and reviews11,15,16,17,18 comparing various algorithms to generate conformers. The metric used throughout is an algorithm’s ability to identify a conformation with low structural RMSD (which is surely important) to an a priori known bio-active conformation but neglecting the algorithm’s ability to energy-rank this conformation. For further information we refer to a very recent and comprehensive review by Hawkins11 on “Conformation Generation: The state of the Art” covering purpose, methods, validation as well as current and future directions in this field. In this paper we deal with step two, the energy prediction. We benchmark force field and semiempirical quantum mechanical methods as well as the speed- and accuracy-optimized


4

Page 5 of 49 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


PBEh-3c density functional19 against PBE020 with D3(BJ) dispersion correction21,22 and the def2TZVP basis set23 with regard to their capability to reliably provide conformer energies, conformer rank-orders and a representative low-energy conformer per molecule on a carefully selected set of 100 pharmaceutical drug substances. We do not compare conformer generator algorithms but selected the BEST algorithm24,25 as implemented in BIOVIA Pipeline Pilot26 based on its speed, ability to also create conformers for non-standard ring systems, and diversity of raw conformers.27 We provide results in settings of increasing order of complexity, starting with gas-phase structures and energies in the first section, then adding solvation free energies for water and finally we go from neutral compounds to the predominant charge states at pH 7.4. We also provide results for the dependence of the energies on the geometry optimization method chosen. We explicitly ignore tautomerism, being aware that the tautomer state may change due to solvation and protonation state.

METHODS, METRICS, DATASET DATASET PREPARATION All calculations were performed on a selection of 100 drug-like molecules from the cytochrome P450 data set that was published along with the XenoSite model by Zaretzki et al.28 using the chemical structures exactly as reported in this paper. The Xenosite dataset consisting of 680 molecules was first filtered down to a molecular weight in the range of 250 g/mol to 500 g/mol which is generally considered to be lead- to drug-like. Then, compounds with more than 13 rotatable bonds were removed. By this unusually high number also some very flexible molecules were kept to challenge the software used in our benchmark. The molecules were grouped according to their numbers of rotatable bonds, and


5


Page 6 of 49

randomly about every fourth molecule from each rotatable bond group was selected yielding the final dataset of 100 molecules shown in Scheme 1. By a recent publication of the Kirchmair group29 we lately became aware of issues with wrong chemical structures in the Xenosite dataset. We thoroughly checked all structures and found deviations with regard to stereochemistry for 15 compounds, namely 5, 9, 13, 25, 28, 31, 38, 39, 41, 51, 82, 84, 86, 93 and 95. Additionally, in compound 2 the para fluorine in the phenyl ring is missing, as is a double bond in the indole ring of 45, and compound 48 carries an alcohol instead of an acid. We have corrected the chemical structures and provide them as structure files in Supporting Information.


6

Page 7 of 49 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

32

33

34

35


7


Page 8 of 49

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

61

62

63

64

65

66

67

68

69

70


8

Page 9 of 49 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


71

72

73

74

75

76

77

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

99

100

78

97

98

Scheme 1 Dataset of 100 molecules.


9


Page 10 of 49

3D coordinates were created in Pipeline Pilot26 after addition of hydrogen atoms using the components “3D Coordinates” and “Add Hydrogens” with standard settings. A maximum of 1000 conformers were generated per molecule by the BEST algorithm24,25 implemented in Pipeline Pilot, applying initially a high absolute energy filter of 50 kcal/mol to circumvent the risk of removing novel distinct conformations due to potentially sub-optimal internal clustering of the conformer algorithm. In a last step, the conformers were filtered by BEST relative energy values keeping up to 100 of the lowest conformers but discarding those exceeding 20 kcal/mol. The BEST algorithm had been chosen from the portfolio of available software at Bayer because it had the best compute cost plus diversity benefit ratio in an internal evaluation27 and especially because it is able to create non-standard and macrocycle ring conformers. CALCULATION METHODS For the first part of our benchmark, molecules were kept in neutral state except for molecule 42 which is a quaternary amine. All conformers were geometryoptimized in gas-phase with the recently published semiempirical method GFN-xTB30. Briefly citing the description of the authors, GFN-xTB “is related to the self-consistent density functional TB scheme and mostly avoids element-pair-specific parameters. The parametrization covers all spd-block elements and the lanthanides up to Z = 86 using reference data at the hybrid density functional theory level. Key features of the Hamiltonian are the use of partially polarized Gaussian-type orbitals, a double-zeta orbital basis for hydrogen, atomic-shell charges, diagonal thirdorder charge fluctuations, coordination number-dependent energy levels, a noncovalent halogen-bond potential, and the well-established D3 dispersion correction.” Gas-phase energy calculations were then performed with the widely used force field MMFF9431 and the more recent force field OPLS332 in the Schrödinger suite 2017-233, the semiempirical methods AM134 as available from the VAMP component35 in Pipeline Pilot and


10

Page 11 of 49 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


DFTB-D336,37,38,21 implemented in the software DFTB+ version 1.2,39 which is another tightbinding approach, and finally with the efficiency-optimized PBEh-3c density functional.19 As benchmark reference we used PBE020 with D3(BJ) dispersion correction21 and the def2-TZVP basis set23. PBEh-3c and PBE0 calculations were performed with TURBOMOLE 7.1.40,41 Geometry optimizations and TURBOMOLE setup were carried out with the software ancopt42,43. Free energies of solvation with solvent water were obtained for the neutral conformers in the respective solvation models of the software as applicable. The coordinates used throughout were optimized with GFN-xTB using the implemented GBSA model.44 Energy calculations for MMFF94 and OPLS3 were done with the GBSA-model of the Schrödinger suite, AM1 solvation was achieved with the SCRF model,45 and PBEh-3c and PBE0 free energies of solvation were calculated with COSMO-RS46,47,48 continuum solvation in COSMOtherm 160149 with the BP86 functional,50 def2-TZVPD basis set and m3 grid. Since there is no solvation model for DFTB+, it was excluded from the analysis. For the third setting of the benchmark, predominant charge states at pH 7.4 (pH in blood plasma) for all molecules were calculated by the pKa module51 of the software ADMET PREDICTOR 7.1 by Simulations Plus.52 All conformers were again optimized with GFN-xTB using the implemented GBSA model and energies and solvation free energy corrections were applied as described in the previous paragraph. For the settings just described we always optimized the coordinates of all conformers with the semi-empirical method GFN-xTB. We were aware that any choice of one method would introduce a certain bias against all other methods. Therefore, in a fourth setting we provide results on the dependence of the energy rankings on geometry optimization. Here, geometry optimizations were additionally carried out with OPLS3 with native solvation model followed by


11


Page 12 of 49

PBE0-D3(BJ)/def2-TZVP energies with COSMO-RS solvation free energy corrections added as well as with PBEh-3c in gas-phase and PBE0-D3(BJ)/def2-TZVP gas-phase energies. There are some method-intrinsic limitations. Due to missing parameters for bromine we cannot provide DFTB+ results for compounds 7 and 25. Likewise calculations for 68 were not possible with DFTB+ and MMFF94 due to missing parametrizations for boron.

Figure 1 Statistical metrics

METRICS For each molecule all conformers were labeled by a counter. In a first step, we calculated the relative energies of each conformer to the lowest energy conformer in each method. Then we compared the respective values of each labeled conformer for the method at hand with the values calculated by the benchmark method, be it the relative energy or the rank.


12

Page 13 of 49 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


By this we derived the metrics per molecule per method by aggregating over all conformers of a molecule. Those metrics were then further aggregated on a method level, i.e. we report the median and maximum values for all over all molecules considered in our analysis in Tables S2 to S4 of Supporting Information (information on molecules omitted from statistics for some reason will be given in the respective sections). The statistical metrics (see Figure 1) can be categorized into •

standard error measures as are mean absolute error (MAE), RMSE and maximum error (MAX),

•

the correlation metrics Pearson correlation, and the Spearman coefficient, which is a metric of rank-order conservation, and

•

DeltaEmin, the relative PBE0 energy of the conformer ranked lowest according to the method under investigation.

In the following, we will consistently focus on median values to avoid bias due to any outliers and provide the fairest comparison. Nevertheless, the box plots presented give also information on the value spread and distribution. It is important to keep in mind that all results compare calculated values with calculated values of the benchmark method. The benchmark method itself has a weighted total mean absolute deviation of 3.6 kcal/mol similar to the other methods tested in a benchmark study by Goerigk et al.54 The relative error for conformer energies for a dispersion-corrected hybrid functional with sufficiently large basis set is known to be considerably lower due to error cancellation. BOXPLOTS The boxplots were generated with Microsoft Excel 2016. The median values shown are inclusive medians. The box stretches from the lower quartile to the upper quartile of the data. In line with the standard settings for boxplots in Microsoft Excel 2016 the whiskers do not


13


Page 14 of 49

represent minimum and maximum values of the dataset, but stretch from the lowest value within a range of 1.5 times the interquartile range (IQR) from the lower quartile to the highest value within a range of 1.5 times the IQR from the upper quartile. STATISTICS We applied pairwise t-Test statistics to the six error metrics using a short R script.53 The assumption of normal distribution is not strictly fulfilled (cf. data distributions in the box plots). P-values are adjusted for the fact of multiple testing of each method against all others. We set the level of significance, alpha, to 0.001 for the probability of the null hypothesis. We test the H0 that two methods have equal error rates or correlation coefficients. Pairwise adjusted pvalues are provided as Supporting Information in Tables S5, S6, S7.

RESULTS In this work we provide benchmark results for relative conformer energies of the MMFF94 and OPLS3 force fields, semiempirical AM1, DFTB+ and GFN-xTB methods and PBEh-3c density functional theory calculations compared to the density functional PBE0D3(BJ)/def2-TZVP that itself has been benchmarked by Grimme et al.54, 19 against a great variety of other QM and DFT methods and proven to be highly reliable. Similar results were obtained in a benchmark study55 on amino acid conformations. There, PBE0-D3(BJ) was on par with other GGA, meta-GGA and hybrid functionals, but weaker than double-hybrids and coupled-cluster approaches. Basis set variability tests in this study clearly show that the triple-zeta level def2TZVP is sufficient. Our benchmark was run at three levels of complexity. At the first complexity level, conformational energies were calculated for neutral molecules (with the exception of 42) in the gas phase. In the second setting, the calculations were repeated in water solvation, which bears the risk of adding noise due to differences in the solvation models. Considering that 47 out of the


14

Page 15 of 49 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


100 molecules in our dataset are charged at physiological pH and the same is true for similar percentages of the compound decks of pharmaceutical companies, the final approach was to calculate conformational energies in water with the molecules at their predominant charge state at pH 7.4 which is the pH in plasma. Charge states were predicted by ADMET predictor52 from Simulations Plus51 before conformer generation. DATASET Scheme 1 shows all structures in their predominant charge state at pH 7.4. Out of the 100 molecules, 32 are positively charged at pH 7.4, 14 are negatively charged, and molecule 10 is a zwitterion. Figure 2 gives an overview of the distribution of numbers of raw conformers. Most of the molecules have more than 10 raw conformations. 59 molecules in gas-phase and 60 solvated charged ones have more than 100 raw conformations that finally were filtered down to the 100 with lowest energy. The overall distributions are very similar for the two cases neutral gas-phase and solvated charged, but for the individual molecules the numbers are sometimes significantly different, as discussed later. The seven molecules 14, 15, 17, 18, 19, 80, 94 with less than five conformers are omitted from statistical analyses in order to achieve conclusive correlations but are discussed in detail in the respective sections.


15


35

Molecules

30 25 20 15 10 5 0 5

10

20

50 100 Bin (Conformers)

200

500

1000

Figure 2 Binned numbers of molecules with a certain number of raw conformers for the settings neutral/gas-phase (grey) and solvated/predominant charge state (dark grey). Bin labels are maximal numbers of conformers in respective bin.

40 35 30 Molecules

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

Page 16 of 49

25 20 15 10 5 0 5

10 15 20 Bin (Gas-phase energy range in kcal/ mol)

>20

Figure 3 Binned numbers of molecules with respective gas-phase energy ranges for the maximally 100 conformers considered. Bin labels are maximal energy in kcal/mol in respective bin. Bars are PBE0 (blue), PBEh-3c (dark blue), AM1 (violet), DFTB+ (light blue), GFN-xTB (orange), MMFF94 (green), OPLS3 (dark red), OPLS3 optimized (grey).


16

Page 17 of 49 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 distribution of gas-phase energy ranges, i.e. the difference between the highest and minimum conformer energies per molecule (Figure 3), is quite different for the different methods, with a slight grouping into DFT (PBE0, PBEh-3c) plus AM1, tight-binding QM and force field methods. Force field energy calculations show a wider spread and a strong peak for > 20 kcal/mol, i.e. they extend even outside the range of the BEST raw conformer energies. After geometry optimization (starting from the GFN-xTB coordinates) the bin profiles of the force fields (as shown exemplary for OPLS3 optimized coordinates as grey bars in Figure 3) are similar to the QM profiles, which hints to suboptimal bond lengths and angle parameters being responsible for the shifts to higher values. We are aware that this in the end might lead to artificially weaker performances of the force field methods. We analyse this effect in more detail in section “Influence of Geometry Optimization”. NEUTRAL MOLECULES, GAS PHASE For our initial benchmark setting, the BEST-generated conformers were optimized with GFN-xTB. Energy calculations were then performed by the other methods with the intention to exclude any effects due to different coordinates and to focus on energy-ranking. PBE0-D3(BJ)/def2-TZVP was used as reference energy method. Obviously, in a practical setting one would not optimize geometries with a QM method and do energy calculations with a force field.


17


A)

B)

Page 18 of 49

C)

Figure 4 Boxplots for gas-phase values for A) mean absolute error MAE, b) root mean square error RMSE, and C) maximal error MAX over 93 molecules, in units of kcal/mol. Bars are BEST (blue), OPLS3 (dark red), MMFF94 (green), AM1 (violet), DFTB+ (light blue), GFNxTB (orange), PBEh-3c (dark blue).

In the following, in the text we provide median values of the aggregated medians per molecule for the performance of each method, additionally to the box plots which also give the variability over the data set. As expected (Figure 4), PBEh-3c shows extremely low MAE, RMSE and MAX errors of 0.67, 0.77 and 1.89 kcal/mol and also error variability. It performs significantly better (see Supporting Information Table S4)

on all metrics compared to almost all methods. This confirms the

benchmark results versus TPSS-D356 with def2-TZVP basis set from the original publication. The other extreme, as expected, is our negative benchmark BEST energies (which are based on BEST coordinates in contrary to GFN-xTB coordinates as for the other methods) where the errors are twice as high as for the second worst method (again statistically significant for the comparison of to AM1, DFTB, GFN-xTB and PBEh-3c). Obviously, those energies should not be used. BEST data are therefore excluded from all further analyses.


18

Page 19 of 49 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


GFN-xTB outperforms all other methods except PBEh-3c in error measures (significant for comparison to AM1 and DFTB), followed by the other tight-binding ansatz DFTB+. AM1 has errors more in the range of the force fields, and OPLS3 is superior to MMFF94.

A)

B)

C)

Figure 5 Boxplots for gas-phase values for A) Pearson cc, B) Spearman cc, and C) DeltaEmin over 93 molecules, in units of kcal/mol. Bars are OPLS3 (dark red), MMFF94 (green), AM1 (violet), DFTB+ (light blue), GFN-xTB (orange), PBEh-3c (dark blue).

The overall correlation of the conformer energies calculated with the method observed to the benchmark conformer energies, expressed by the Pearson correlation coefficient, as well as correct rank-ordering of the conformers (even if absolute energies are not described well) described by the Spearman rank coefficient are important metrics for the identification of relevant representative conformers, e.g. in the context of pharmacophore modeling or ligand strain energies. In both disciplines (see Figure 5) AM1 is extremely weak with about random Pearson and Spearman cc’s of 0.51 and 0.44. DFTB+ and OPLS3 are on par, but clearly outperformed by GFN-xTB with 0.86 and 0.79, respectively. PBEh-3c shows almost perfect correlation.


19


Page 20 of 49

Generally we find that for conformer sets where the energy range exceeds about 12-15 kcal/mol the Pearson correlation coefficient is at least 0.7. Smaller energy ranges on the other hand provide no indication on the correlations to be expected. Finally, DeltaEmin, the relative PBE0 energy of the conformer ranked first by the method investigated, is a measure of how far off on median the lowest energy conformer is. GFN-xTB and OPLS3 are clearly below 1 kcal/mol (0.39 and 0.54; difference not statistically significant), DFTB+ and MMFF94 with 1.24 and 1.2 show higher errors and especially AM1 with 2.36 kcal/mol is significantly off. In the following we will focus on OPLS3 as the better performing force field and GFN-xTB as the best-performing semiempirical QM approach with relatively low computing costs (no significant difference in performance), and error-wise we will focus on the MAX error as a worst-case scenario for complicated molecules and DeltaEmin as a measure of how reliable our minimum energy conformer selected by a certain method is. Figure 6A provides a binning of the numbers of molecules with a certain MAX error and Figure 6B a binning for the DeltaEmin, each time in kcal/mol. Plots for the other methods are given in Figure S4 of Supporting Information.


20

40 35 30 25 20 15 10 5 0

Molecules

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


Molecules

Page 21 of 49

2

A)

4 6 8 10 >10 Bin (MAX in kcal/mol)

B)

70 60 50 40 30 20 10 0 1 2 3 4 5 >5 Bin (DeltaEmin in kcal/mol)

Figure 6 Binning of A) maximum error distribution MAX and B) energy deviation DeltaEmin for neutral molecules in the gas phase, GFN-xTB in grey and OPLS3 in dark grey. Binned are numbers of molecules and bin labels are maximal errors in kcal/mol in respective bin.

Whereas GFN-xTB shows a Gaussian-like distribution with a maximum at 6 kcal/mol, OPLS3 has a strong peak at > 10 kcal/mol, i.e. a higher risk to have a wrong ranking of conformers. On the other hand, the binnings for DeltaEmin, the energy error for the best-ranked conformer by each method, are very comparable, with 83 (GFN-xTB) and 78 (OPLS3) molecules out of 93 better than 3 kcal/mol. Taking all measures into account, GFN-xTB is the winner in gas-phase (except for PBEh-3c). It provides the lowest errors, the highest correlation coefficients and the lowest risk of being far off. The difference to OPLS3 nevertheless is not statistically significant. As said, the seven molecules 14, 15, 17, 18, 19, 80, 94 with less than five conformers were omitted in all statistics. Nevertheless, we kept them in the dataset since they allow for detailed inspection of the methods’ capabilities. Molecule 94 has two sterically induced energyequivalent conformations (by the way, most conformer generators are not able to create both from one initial coordinate set) as calculated by all methods, the largest deviation being 0.02 kcal/mol for OPLS3.


21


Page 22 of 49

The only difference between the two conformers of 14 are the orientations of the methyl group and the hydrogen which are almost equivalent in energy with a delta of less than 0.2 kcal/mol, except for AM1 with 0.7 kcal/mol for the PBE0 preferred conformer. The outlier is OPLS3 with the methyl trans to the annelated ring being 1.9 kcal/mol higher in energy. Compound 17 has 3 conformations, 2 of them enantiomeric. OPLS3 and GFN-xTB prefer the singular conformation, whereas all other methods prefer the enantiomeric pair, but with a difference of less than 1 kcal/mol. Two mirror image conformers are generated for 18 which are energy-equivalent in all methods. 19, due to its additional chlorine atom should yield two pairs of mirror images, but we find only three conformers. All methods see the pair equivalent, but AM1 and DFTB+ provide the wrong ranking. GFN-xTB under- and MMFF94 overestimates the difference by 1.5 and 1.6 kcal/mol. Compound 80 has 2 mirror image pairs separated by 1.8 kcal/mol (PBE0) and overestimated by GFN-xTB and MMFF94 by 1.6 and 4.1 kcal/mol. AM1 again provides the wrong rank ordering. Finally, for 15 we find only 1 conformer. The methyl group is kept inplane (even if starting minimizations from 90° dihedral angle) and all rings are flat. With molecule 27 the dataset contains a macrocycle which we expected to be markedly challenging. Nevertheless, all methods succeeded in differentiating between low and high energy conformers with Spearman cc’s of 0.6 or higher. The best result apart from PBEh-3c is obtained by GFN-xTB with a Spearman cc of 0.81. Although GFN-xTB and OPLS3 perform well in general, some molecules appear to be challenging for them. The Spearman cc is below 0.5 for 11 molecules when using GFN-xTB and for 14 molecules when using OPLS3 in gas phase but for different molecules except for 44 and 99. We do not see any structural patterns other than a slight tendency of GFN-xTB to be weaker for molecules with more than one methoxy group (e.g. outlier 9 with Spearman of -0.22 contains


22

Page 23 of 49 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


four methoxy groups). Multiple methoxy groups result in multiple about equi-energetic conformers which are difficult to rank. The other structural feature about molecule 9 is its mix of aliphatic and aromatic bonds in the central ring. Even the two phosphorus containing compounds behave differently with 51 having a high and 97 having a weak Spearman cc for GFN-xTB. The boronic acid 68 is ranked weakly by OPLS3 (no parameters for MMFF94). OPLS3 (and equally MMFF94) underperform for the highly symmetrical 3 and the zwitterion 10 in its uncharged state. Details on all Spearman ranks for all compounds and methods can be found in Figure S1 of Supporting Information. Finally, we report on an occasional problem whenever switching between the QM world which is solely based on atom coordinates and the force field world requiring connection tables including bond definitions. During geometry optimization, bond lengths and angles can leave the typical ranges e.g. via conjugation effects or via intramolecular hydrogen bonds that lead to formal tautomerism. In the beginning, we used a Pipeline Pilot node based on a Maximum Weighted Matching algorithm for nonbipartite graphs57 to reassign bond orders to the 3D coordinates after the GFN-xTB optimization. In a couple of cases, however, the bond orders were assigned incorrectly leading to faulty molecular structures and charges. We thus decided to consider the bond orders to be unchanged throughout the calculations and mapped back the input bond orders onto the coordinates. This in a low number of cases neglects actual changes in bond orders and caused unrealistically high force field energy values. After all, however, considering the bond orders as stable reduced the amount of structural errors to a minimum and those high energies fall out of the median value statistics we provide. Though we are aware of this issue, we accept it for the sake of an automated process.


23


Page 24 of 49

NEUTRAL MOLECULES, WATER SOLVATION In this setting, the BEST-generated conformers were optimized with GFN-xTB in solvent water and energy calculations were performed with the other methods always using the solvent model of the method. PBEh-3c and PBE0-D3(BJ)/def2-TZVP are solvation free energy corrected by COSMO-RS, the latter was used as reference. DFTB+ had to be skipped due to missing solvent model. Most obviously (see Figure 7), all error values rise for GFN-xTB, i.e. for MAE from 1.71 to 2.23, for RMSE from 2.13 to 2.80 and for MAX from 5.21 to 7.14, and also for AM1, whereas the errors for the force fields decrease both by at least 0.3 for all measures. The distributions of the medians as shown in the box plots become more similar for all methods compared to the gasphase distributions. All data presented in the following suggest that this is due to the fact that GFN-xTB and AM1 were parameterized on gas-phase data. MMFF94 and OPLS3Error! Bookmark not defined. were parametrized on experimental and calculated gas-phase data with regard to geometries, but apply fitting schemes for electrostatics to reproduce charge distributions, dipole moments and energies of the liquid state. The publication on OPLS3 (which is an acronym for “Optimized Potentials for Liquid Simulations”) explicitly states that “torsional terms fit to the gas phase potential surface can be transferred to the condensed phase with minimal errors. This hypothesis requires only that the change in potential energy as a function of torsion angle is well represented by the nonbonded interaction terms in the force field.” I.e., stretch, bend and torsional parameters are derived from gas-phase QM but the non-bonded terms for electrostatics and van der Waals-parameters are derived from gas-phase and condensed phase experimental data.58 Atom-centered partial CMA1-BCC charges59 are applied which are scaled by bond charge correction terms60 to minimize errors to experimental absolute solvation free energies. In a similar way MMFF94 was parameterized based on QM-derived bond lengths,


24

Page 25 of 49 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


angles and torsions but with non-bonded interactions scaled in a way to mimic “the enhancement of charge distribution due to molecular polarizability”31 in condensed phase by so-called effect pair potentials.61

A)

B)

C)

D)

E)

F)

Figure 7 Boxplots over 93 molecules in water solvation, in units of kcal/mol for the error plots; A) mean absolute error MAE, b) root mean square error RMSE, C) maximal error MAX, D) Pearson cc, E) Spearman cc, and F) DeltaEmin. Bars are OPLS3 (dark red), MMFF94 (green), AM1 (violet), GFN-xTB (orange), PBEh-3c (dark blue). The error values for PBEh-3c MAE (0.67 to 0.69) and RMSE (0.79 to 0.80) stay constant but MAX is reduced from 2.21 to 1.89, which is caused by the different coordinates compared to gas-phase. Note, that we apply identical solvation free energy corrections for PBEh-3c and PBE0.


25


Page 26 of 49

Overall, based on error measures GFN-xTB is the best-performing method (but only significantly better compared to AM1 on the 0.01% level; see Supporting Information Table S6), even though the force fields catch up. This can also be seen from the DeltaEmin and MAX error dstributions for all methods provided in Supplementary Information Figure S5. The picture changes looking at the correlation coefficients, where both force field methods show unchanged high correlations, whereas the values for GFN-xTB drastically decrease from 0.86 to 0.51 for Pearson and even more from 0.79 to 0.45 for Spearman cc (significant for GFN-xTB vs. OPLS3). AM1 stays as weak as in gas-phase. OPLS3 has 14 compounds with Spearman cc < 0.5 as in gas-phase, seven of which are identical and seven new. The drastic decline of rank-ordering capability of GFN-xTB is reflected in now 41 compounds with Spearman cc < 0.5, six of them identical to the overall eight from gas-phase (detailed data for all compounds and methods can be found in Supplementary Information Figure S2). Overall, this indicates that a relatively modest increase in error measures drastically changes the orders of the conformers, which indicates that the errors might be more systematic for force fields than for the semiempirical methods. Finally, the pronounced change in DeltaEmin also indicates that the rank-ordering of the conformers is highly influenced by going from gas-phase to solvation. Overall, one can assume that the GFN-xTB solvent model is the main source of the deterioration seen. We therefore compared the median absolute differences of the free energies of solvation from GFN-xTB and OPLS3 to COSMO-RS. The median values for the 93 minimum energy conformers are 5.5 and 3.3 kcal/mol for GFN-xTB and OPLS3, whereas the median values over all 7709 conformers are 13.7 and 3.5 kcal/mol. The correlation coefficients to COSMO-RS for the free energies of


26

Page 27 of 49 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


solvation for the 93 mimimum energy conformers are 0.74 and 0.77 for GFN-xTB and OPLS3, and for all conformers are 0.73 and 0.80, respectively. Correlation plots can be found in Supporting Information Figure S7. This solvation free energy deviations between the solvation models of the two methods and COSMO-RS explains the differences in rank-ordering capability.

MOLECULES IN PREDOMINANT CHARGE STATE, WITH WATER SOLVATION In the final setting, we used ADMET predictor52 to identify the predominant charge state for the molecules at pH 7.4 (plasma pH). Out of the 100 molecules, 48 are now non-neutral. Charge state prediction was performed for the tautomeric form as shown in Scheme 1. Conformers were generated with the BEST algorithm after protonation, followed by optimization with GFN-xTB with its inherent solvent model and energy calculations with the other methods, always using the respective solvent model. PBE0-D3(BJ)/def2-TZVP with COSMO-RS solvation free energy correction was used as reference energy method. DFTB+ had to be skipped due to missing solvent model. For about one third of the charged structures, we observe relevant differences of the numbers of raw conformers per molecule (i.e. the numbers before filtering down to maximally 100) compared to the neutral state, strongly connected to specific chemical functionalities. Whereas carboxylic groups induce slight changes (-20 to 20) of conformer numbers, central aliphatic ring nitrogens and especially piperidine nitrogens raise the numbers by up to 490 additional conformers, the exception being compound 43 (-75 raw conformers). Noncyclic amines change the numbers more moderately and in both directions. This can be rationalized by the fact that the acidic groups and non-cyclic amines are terminal yielding a low number of alternate


27


Page 28 of 49

conformations. On the other hand, the positively charged nitrogen-containing rings are central, acting as lever arms for muultiple diverse conformations and additionally result in enantiomers adding many new conformers to the set.

A)

B)

C)

D)

E)

F)

Figure 8 Boxplots over 93 molecules in their predominant charge state in water solvation, in units of kcal/mol; A) mean absolute error MAE, b) root mean square error RMSE, C) maximal error MAX, D) Pearson cc, E) Spearman cc, and F) DeltaEmin. Bars are OPLS3 (dark red), MMFF94 (green), AM1 (violet), GFN-xTB (orange), PBEh-3c (dark blue).

Nevertheless, going from neutral to predominantly charged conformers as shown in Figure 8 does not change the overall statements made for the comparison of gas-phase versus solvation


28

Page 29 of 49 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


phase value changes, i.e., the median ability to calculate energies and to rank molecules over a larger set is not charge-state dependent. All PBEh-3c values again stay almost constant. All error measures for all methods slightly increase by about 0.1 to 0.4 kcal/mol except for the high increase in OPLS3 MAX error of 1.15 kcal/mol. The errors are always lower for GFN-xTB than for OPLS3 which performs better than AM1 and MMFF94 except for the MAX error where AM1 is better with 8.51 compared to 8.97 kcal/mol of OPLS3, but with a much larger spread of values. Supporting information Table S7 provides information on the significance of the differences described. The changes in correlation coefficients are even less pronounced except for AM1 which goes from weak to bad. Both Pearson and Spearman decrease by 1 to 5 percent except for a slight increase of the GFN-xTB Spearman value. The same applies for DeltaEmin which increases by 0.68, 0.42 and 0.87 kcal/mol for OPLS3, MMFF94 and AM1 but decreases by 0.36 kcal/mol for GFN-xTB. Overall, the differences are not statistically significant except for the pair OPLS3 / AM1. The median Spearman rank for GFN-xTB goes down dramatically from 0.79 for gas-phase to 0.62 for charged molecules in solvent. The reason is a massive increase of molecules with Spearman rank below 0.5 from 11 (gas) to 41 (solv, neutral) and 36 (solv, charged). The corresponding numbers for OPLS3 are 14, 14, and 22 molecules. The detailed analysis of MAX error in Figure 9A and DeltaEmin in Figure 9B (Plots for the other methods are given in Supporting Information Figure S6) if compared to Figure 6 shows a strong shift to both higher MAX errors and lower numbers of molecules with small DeltaEmin for GFN-xTB and a certain shift for OPLS3 compared to gas-phase. GFN-xTB is still better in


29


the former but now worse in the latter measure. Nevertheless, DeltaEmin is below 3 kcal/mol for 71 and 68 molecules out of the 93 considered for GFN-xTB and OPLS3, respectively.

40 35 30 25 20 15 10 5 0

Molecules

Molecules

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

Page 30 of 49

2

A)

4 6 8 10 >10 Bin (MAX in kcal/mol)

B)

70 60 50 40 30 20 10 0 1 2 3 4 5 >5 Bin (DeltaEmin in kcal/mol)

Figure 9 Binning of A) maximum error distribution MAX and B) energy deviation DeltaEmin for solvated molecules in predominant charge state, GFN-xTB in grey and OPLS3 in dark grey. Bin heights are numbers of molecules and bin labels are maximal energy errors in kcal/mol in respective bin.

Quite positively, all methods perform similarly for neutral (previous section) and charged species in solvation. The poor Spearman ranks are about evenly distributed over the two classes. The only zwitterion, namely molecule 10, is an indicator of the strengths and weaknesses of quantum mechanics versus force field. In gas-phase, OPLS3 does not and GFN-xTB almost perfectly ranks the 100 (145 raw) conformers (Spearman: -0.09 vs. 0.86), whereas the ranking of the non-zwitterionic solvated form reflects the general observations for the change in medium, in that OPLS3 Spearman increases to still weak 0.45 and the other decreases to 0.25. For the zwitterionic form, the situation changes again with Spearman cc’s of -0.22 and 0.52. This clearly demonstrates two effects being mixed in the outcome: (i) only a wave function method is able to


30

Page 31 of 49 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


describe the preferred electronic states in gas-phase and polar solvent, whereas (ii) the description of solvent effects for OPLS3 is superior to the GBSA model of the QM method.

INFLUENCE OF GEOMETRY OPTIMIZATION In the first three sections we always provided the same input coordinates as derived by GFN-xTB optimization in gas-phase or with its inherent solvent model to all methods. We always compared single point energies to remove any bias from geometry optimization which would lead to different minima and make energies not comparable anymore. The assumption that is used in the plethora of quantum chemistry publications is that potential energy hypersurfaces of different methods are parallel. In an earlier publication we were able to provide some evidence that this might not be the case12 especially for larger and more globular molecules. Here we test the parallel hypersurface assumption by two scenarios. In the first scenario, we compare the correlation of gas-phase energies for PBEh-3c to PBE0 using either GFN-xTB or PBEh-3c optimized coordinates. As shown earlier, PBEh-3c is the by far best performing method in our benchmark but computationally out of reach for most applications in this area of research, where thousands of molecules have to be processed in short time (PBEh-3c geometry optimizations should be avoided therefore if possible). Consistently, all errors increase and correlations decrease (see Table 1, top) when going from a comparison of PBEh-3c with PBE0 both optimized with GFN-xTB to a comparison of PBEh-3c energies from PBEh-3c optimization with GFN-xTB optimized benchmark energies. But then, doing the comparison with both times PBEh-3c optimized coordinates the errors decrease and correlations increase compared to the purely GFN-xTB coordinate values, indicating that the potential energy surfaces (PES) for PBEh-3c and PBE0 are mire similar than for the GFN-xTB


31


Page 32 of 49

and PBE0. The RMSE, for instance, goes from 0.79 via 1.61 for the “mixed” coordinates to 0.62 kcal/mol. Though all metrics are slightly better for PBEh-3c coordinates, overall PBEh-3c//GFNxTB yields highly reliable conformer geometries and energies and is computationally cheaper by a factor of 100 to 1000. Table 1 Influence of geometry optimization method on statistical metrics. Column one describes how the structure for the method considered is optimized and column two the optimization for the benchmark method PBE0-D3(BJ)/def2-TZVP. The two methods considered are PBEh-3c and OPLS3. Method

Benchmark

MAE

RMSE

Spearman

Pearson

MAX

Delta Emin

Geometry optimization by Energies for PBEh-3c vs. Benchmark GFN-xTB

GFN-xTB

0.67

0.79

0.97

0.98

1.89

0.00

PBEh-3c

GFN-xTB

1.12

1.61

0.87

0.88

5.13

0.18

PBEh-3c

PBEh-3c

0.53

0.62

0.97

0.98

1.52

0.00

Energies for OPLS3 vs. Benchmark GFN-xTB

GFN-xTB

2.4

3.16

0.71

0.7

8.97

1.29

OPLS3

GFN-xTB

2.11

2.95

0.63

0.62

8.19

1.90

OPLS3

OPLS3

1.66

2.00

0.76

0.81

4.97

0.96


32

Page 33 of 49 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


In the second scenario, we compare the solvent-phase predominant charge state calculations for OPLS3 in the same way. Here, the effect of the non-parallel hypersurfaces is more pronounced, which was to be expected due to the closer relationship between the three methods GFN-xTB, PBEh-3c and PBE0. Here, all error measures go down going from GFN-xTB coordinates to OPLS3 coordinates, independent of the origin of the benchmark coordinates. On the other hand, Pearson and Spearman coefficients go down and DeltaEmin up for the “mixed” comparison. But for the pure OPLS3 setting, all errors and DeltaEmin are significantly lower than for the setting with GFN-xTB coordinates, and the cc’s are higher.

COMPUTATIONAL RESOURCES The timings for the different job types given in the following are mean values per conformer derived from the 18 molecules we had to re-calculate due to the mentioned errors in chemical structures. The numbers are typical values with low variances for energy calculations and larger variations for geometry optimizations depending on the starting coordinates due to variability in the numbers of optimization cycles needed (between about 20 to more than 200). Calculations were performed on our compute cluster using three diffenent compute node types applying Intel CPUs with 2.4 to 3.2 GHz and between 4 and 12 cores. Whereas AM1 (Semiempirical QM descriptors) and BEST (Conformation Generator) calculations were run using native Pipeline Pilot components, all other calculations were set up and run via Pipeline Pilot workflows using the component “Program (run on Server)” which starts the actual executables for calculation. The Pipeline Pilot process overhead is nevertheless below 1 % for the actual job submission as reported by the functionality “Show Process Time”. All values are rounded mean values for the calculation of one conformer but based on the values for 1919


33


Page 34 of 49

neutral conformers, 1919 solvated neutral conformers and 2148 solvated charged conformers. Timings are for 1 CPU core, always including Pipeline Pilot processing overhead. GFN-xTB is generally running as multi-thread process. BEST raw conformer generation needed 0.4 s per conformer. Energy calculations for OPLS3 and MMFF94 run for 7.2 s and 1.3 s, respectively, whereas AM1 and DFTB+ run for 0.6 s and 8 s. GFN-xTB needs about 7.7 s. PBEh-3c and PBE0 are by orders of magnitude slower, but timing significantly depends on molecule size, i.e. number of basis functions. Whereas the computing time for the single points is quite similar between OPLS3 and GFNxTB, the timings for geometry optimization differ significantly, with 8.3 s for OPLS3 and 127 s for GFN-xTB. Taking together the overall method performances and the requested compute times, it is advisable to optimize the BEST raw conformers in the predominant charge state and in solvent with OPLS3 as a robust and relatively fast method. The whole process of conformer generation, pKa prediction and energy optimizations of maximally 100 conformers can be performed in about 6 min of CPU time per ligand when using OPLS3 which was tested on a set of roughly 30,000 recent project compounds from the Bayer compound deck.

SUMMARY AND CONCLUSIONS In this work we present a benchmark study on the performance of force field, semiempirical and density functional methods for the calculation of relative conformer energies for 100 compounds of lead-like and drug-like size and properties. Our motivation is to understand the strengths and weaknesses of the different methods beyond the “known knowns” that force fields might fail on


34

Page 35 of 49 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


certain substructures due to missing parameters and may be weak for non-classical interactions or that only a first principles ansatz might be able to correctly predict energies. Our motivation is in particular to identify a method able to reliably describe the set of low-energy conformers of any typical molecule. Coming from pharmaceutical industry, we have the additional constraint that the method has to be fast enough to e.g. calculate between one thousand and ten thousand molecules in maximally half a day on a cluster of one hundred cores, which excludes standard ab initio and density functional methods with large basis sets of at least triple-zeta quality required to give quantitative results. Sets of low-energy conformers are relevant for many applications in pharmaceutical research, from pharmacophore alignments and searches over strain energy prediction for the ligands in protein binding, regioselectivity in synthesis planning, atropisomerism to finally ADMET properties which depend on molecular flexibility, such as membrane permeability or solubility. This publication does not deal with the very first step of the process, namely the generation of the raw conformers. We partially relied on publications by others and additionally did some benchmarks with the software we had at hand. Results will be published elsewhere. We have set up the benchmark in three levels of complexity to be able to isolate the different components influencing the results as much as possible. The biggest influence on making data incomparable are different coordinate sets derived from optimizing the geometries with different methods. Therefore, we always used GFN-xTB optimized coordinates and calculated energies on top of those with all methods. In the first setting we looked into gas-phase energies for neutral molecules (with the exception of 42, a quaternary amine). In the second setting we added water solvation which represents the


35


Page 36 of 49

state we are really interested in. For each method we applied the “native” solvation model of the software used if available, and COSMO-RS for the two density functionals. In the third setting we additionally went from neutral to predominantly charged molecules with water solvation. Finally we looked into the effects of geometry optimization for OPLS3 and PBEh-3c. We explicitly excluded tautomerism from our study since there is no cost-effective and conclusive way to identify the correct tautomer in gas-phase or solvation in the correct charge state. As expected, the “low-cost” density functional PBEh-3c is the clear winner but still too expensive for many of our daily applications. It is an extremely reliable and fast DFT approach with almost perfect correlation to our benchmark PBE0-D3(BJ)//def2-TZVP, which itself was previously shown to be extremely reliable. This is true for molecules in neutral and predominantly charged state. We mention here explicitly, that the effect of solvation cannot be judged since both PBEh-3c and benchmark method have the same solvent free energy corrections. PBEh-3c is strongly recommended by us as the method of choice whenever calculation time is no issue. The recently published semiempirical tight-binding approach GFN-xTB is by far the most reliable of the faster methods in gas-phase, which is exemplified by all error metrics, by overall correlation to benchmark as well as by rank-ordering capability. The second semiempirical method, the widely used NDDO method AM1 should obviously not be used. AM1’s error metrics are slightly superior to the force fields in gas-phase, but worse in condensed phase. The correlation coefficients are always very low and DeltaEmin as the expected error for the lowest AM1-predicted conformer is significantly higher than for all methods. Finally, DFTB+, a tightbinding method, is comparable to OPLS3 in gas-phase but not applicable for typical problems


36

Page 37 of 49 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


due to the missing solvent model and missing atom pair parameterizations for boron and bromine. Between the two force fields there is a clear differentiation in performance. Nevertheless, the quite old MMFF94 is better than expected and overall outperforms AM1. OPLS3 in our typical application scenario as is water solvation and predominant charge state shows superior Pearson correlation and rank-ordering capability as expressed by median Spearman coefficient and numbers of molecules with Spearman lower 0.5, combined with lower DeltaEmin, compared to the semiempirical method GFN-xTB. In line with our expectations when we started this study, GFN-xTB is superior to all other fast methods in gas-phase. Somewhat surprisingly, OPLS3 outperforms GFN-xTB in the condensed state. Both methods apply GBSA solvent models which depend on atomic radii and charges. GFN-xTB is parameterized on computed gas-phase data, whereas non-bonded terms and especially the atomic charges of OPLS3 and MMFF94 are fitted against experimental solvation free energies to minimize errors in condensed state. Probably for this reason the force fields developed for condensed phase perform weaker in gas-phase than condensed phase. However, this does not explain why OPLS3 and COSMO-RS solvation free energies are stronger correlated and show significantly lower median errors than the GFN-xTB solvation free energies. Further studies are necessary but out of scope for this publication. Overall, OPLS3 is the obvious choice for many applications in molecular modeling, since it combines short calculation times with high quality conformations. We had expected that the force fields would cause strong outliers due to not appropriately parameterized chemical features, but in the 100 compounds of our study, we saw no particular substructure patterns or


37


Page 38 of 49

functional groups where OPLS3 failed. In the case of such outliers, at least if caused by wrong torsional barriers, the Schrödinger software suite provides a process, the so-called Force Field Builder, that allows to fit parameters for torsions unknown to the force field. We plan to test how such specific fitting influences the statistics for the current dataset as well as for outliers due to poor parametrization.


38

Page 39 of 49 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


AUTHOR INFORMATION Corresponding Author * E-mail: [email protected] Present Addresses $ Master student, University of Hamburg, Germany. # Application manager, University of Wuppertal, Germany. Author Contributions The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Funding Sources There was no funding involved for the performance of this project. ACKNOWLEDGMENT We thank Stefan Grimme for access to the ancopt and GFN-xTB software and also many for fruitful discussions. We especially thank the reviewers for many requests and comments that significantly enhanced this publication.

ASSOCIATED CONTENT Supporting Information. The following files are available free of charge. Additional tables with details to the molecules and comprehensive metrics and statistical data, as


39


Page 40 of 49

well as additional charts. (PDF) Smiles strings and identifiers of the dataset: set100.2D.csv (CSV) Zipped Structure files (SDF) with 3D coordinates: a) neutral raw conformers from BEST algorithm: BEST_RawConformers_neutral.sdf (SDF) b) charged raw conformers from BEST algorithm: BEST_RawConformers_charged.sdf (SDF c) GFN-xTB optimized conformers (gas-phase, neutral): GFNXTB_gas_neutral.sdf (SDF) d) c) GFN-xTB optimized conformers (water, charged): GFNXTB_solv_charged.sdf (SDF)


40

Page 41 of 49 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 Of Contents Graphic:

REFERENCES

1

Alex, A.; Millan, D. S.; Perez, M.; Wakenhut, F.; Whitlock, G. A. Intramolecular Hydrogen

Bonding to Improve Membrane Permeability and Absorption in Beyond Rule of Five Chemical Space Med. Chem. Commun., 2011, 2, 669-674. 2

Over, B.; McCarren, P.; Artursson, P.; Foley, M.; Giordanetto, F.; Grünberg, G.; Hilgendorf,

C.; Lee, M. D.; Matsson, P.; Muncipinto, G.; Pellisson, M.; Perry, M. W. D.; Svensson, R.; Duvall, J. R.; Kihlberg, J. Impact of Stereospecific Intramolecular Hydrogen Bonding on Cell Permeability and Physicochemical Properties J. Med. Chem., 2014, 57, 2746-2754.


41


3

Page 42 of 49

Doak B.C.; Over B.; Giordanetto F.; Kihlberg J. Oral Druggable Space Beyond the Rule of 5:

Insights from Drugs and Clinical Candidates Chem. Biol. 2014, 21, 1115-1142. 4 Leung, S. S. F.; Mijalkovic, J.; Borrelli, K. Jacobson, M. P. Testing Physical Models of Passive Membrane Permeation J. Chem. Inf. Model., 2012, 52, 1621-1636. 5

Menichetti, R.; Kanekal, K. H.; Kremer, K.; Bereau, T. In Silico Screening of Drug-

Membrane Thermodynamics Reveals Linear Relations Between Bulk Partitioning and the Potential of Mean Force, J. Chem. Phys. 2017, 147, 125101; arXiv:1706.02616. 6

Finkelmann, A. R.; Göller, A. H.; Schneider, G. Robust Molecular Representations for

Modeling and Design Derived from Atomic Partial Charges Chem. Commun., 2016, 52, 681-684. 7 Finkelmann, A. R.; Göller, A. H.; Schneider, G. Site of Metabolism Prediction Based on ab initio Derived Atom Representations ChemMedChem, 2017, 12, 606-612. 8 Fu, Z.; Li, X.; Miao, Y.;. Merz, K. M. Conformational Analysis and Parallel QM/MM X-ray Refinement of Protein Bound Anti-Alzheimer Drug Donepezil J. Chem. Theory Comput., 2013, 9, 1686-1693. 9

Sitzmann, M.; Weidlich, I. E.; Filippov, I. V.; Liao, C.; Peach, M. L.; Ihlenfeldt, W.-D.;

Karki, R. G.; Borodina, Y. V.; Cachau, R. E.; Nicklaus, M. C. PDB Ligand Conformational Energies Calculated Quantum-Mechanically J. Chem. Inf. Model., 2012, 52, 739-756. 10

Perola, E.; Charifson, P. S. Conformational Analysis of Drug-Like Molecules Bound to

Proteins: An Extensive Study of Ligand Reorganization upon Binding J. Med. Chem., 2004, 47, 2499-2510.


42

Page 43 of 49 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


11 Hawkins, P. C. D. Conformation Generation: The State of the Art J. Chem. Inf. Model., 2017, 57, 1747-1756. 12 Ehrlich, S.; Göller, A. H.; Grimme, S. Towards full Quantum-Mechanics-based Protein– Ligand Binding Affinities ChemPhysChem, 2017, 18, 898-905. 13 Klebe G. Applying Thermodynamic Profiling in Lead Finding and Optimization Nature Rev. Drug Disc., 2015, 14, 95–110. 14 Lovering, F.; Bikker, J.; Humblet, C. Escape from Flatland: Increasing Saturation as an Approach to Improving Clinical Success J. Med. Chem., 2009, 52, 6752–6756. 15 Kirchmair, J.; Laggner, C.; Wolber, G.; Langer, T. Comparative Analysis of Protein-Bound Ligand Conformations with Respect to Catalyst's Conformational Space Subsampling Algorithms J. Chem. Inf. Model., 2005, 45, 422-430. 16 Friedrich, N.-O.; Meyder, A.; de Bruyn Kops, C.; Sommer, K.; Flachsenberg, F.; Rarey, M.; Kirchmair, J. High-Quality Dataset of Protein-Bound Ligand Conformations and Its Application to Benchmarking Conformer Ensemble Generators J. Chem. Inf. Model., 2017, 57, 529-539. 17 Schwab, C. H. Conformations and 3D Pharmacophore Searching Drug Discov. Today: Technol., 2010, 7, e245 - e253. 18 Agrafiotis, D. K.; Gibbs , A. C.; Zhu , F.; Izrailev, S.; Martin, E. Conformational Sampling of Bioactive Molecules: A Comparative Study J. Chem. Inf. Model., 2007, 47, 1067-1086.


43


Page 44 of 49

19 Grimme, S.; Brandenburg, J. G.; Bannwarth, C.; Hansen, A. Consistent Structures and Interactions by Density Functional Theory with Small Atomic Orbital Basis Sets J. Chem. Phys., 2015, 143, 054107. 20 Adamo, C.; Barone, V. Toward Reliable Density Functional Methods Without Adjustable Parameters: The PBE0 Model J. Chem. Phys. 1999, 110, 6158– 6170. 21 Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate ab initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 elements H-Pu J. Chem. Phys., 2010, 132, 154104. 22 Grimme, S.; Ehrlich, S. Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory J. Comput. Chem., 2011, 32, 1456-1465. 23 Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. RI-MP2: Optimized Auxiliary Basis Sets and Demonstration of Efficiency Chem. Phys. Lett. 1998, 294, 143-152. 24 Smellie, A.; Kahn, S. D.; Teig, S. L. Analysis of Conformational Coverage. 1. Validation and Estimation of Coverage J. Chem. Inf. Comput. Sci., 1995, 35, 285-294. 25 Smellie, A.; Kahn, S. D.; Teig, S. L. Analysis of Conformational Coverage. 2. Applications of Conformational Models. Chem. Inf. Comput. Sci., 1995, 35, 295-304. 26 Pipeline Pilot version 16.5.0.143, Server version 17.1.0.115, Dassault Systemes Biovia Corp. 2016.


44

Page 45 of 49 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


27

Data not shown: we compared ROTATE (Molecular Networks, v. 1.62), ROTATE combined

with CORINA (Molecular Networks, v 3.491) for inclusion of ring conformers, CRESSET conformer algorithm (Cresset, executable xedex, release 2011), MACROMODEL (v103012, maestro 9.7, 2014-1) with different algorithm settings (Schrödinger), FAST, BEST, CAESAR, Systematic Search, Random Search, Boltzmann Jump (Biovia, Pipeline Pilot v.8.5). Diversity was assessed by measures like conformer RMSD distributions or molecular volume diversity. From the two methods with highest structural diversity, macromodel was slower by factors of 10 to 1000 depending on the molecule under consideration. 28

Zaretzki, J.; Matlock, M.; Swamidass, S. J. XenoSite: Accurately Predicting CYP-Mediated

Sites of Metabolism with Neural Networks J. Chem. Inf. Model., 2013, 53, 3373–3383. 29

De Bruyn Kops, C.; Friedrich, N.-O., Kirchmair, J. Alignment-Based Prediction of Sites of

Metabolism, J. Chem. Inf. Model., 2017, 57, 1258-1264. 30 Grimme, S.; Bannwarth, C.; Shushkov, P. A Robust and Accurate Tight-Binding Quantum Chemical Method for Structures, Vibrational Frequencies, and Noncovalent Interactions of Large Molecular Systems Parametrized for All spd-Block Elements (Z = 1 - 86) J. Chem. Theory Comput., 2017, 13, 1989-2009. 31

Halgren T.A. Merck Molecular Force Field. I. Basis, Form, Scope, Parameterization, and

Performance of MMFF94 J. Comput. Chem. 1996, 17, 490-519.


45


Page 46 of 49

32 Shivakumar, D.; Harder, E.; Damm, W.; Friesner, R. A.; Sherman, W. Improving the Prediction of Absolute Solvation Free Energies Using the Next Generation OPLS Force Field J. Chem. Theory Comput., 2012, 8, 2553-2558. 33 Schrödinger Suite 2017-02, Schödinger LLC 2017. 34 Dewar, M.J.S.; Zoebisch, E.G.; Healy, E.F.; Stewart, J.J.P. AM1: A New General Purpose Quantum Mechanical Molecular Model J. Am. Chem. Soc. 1995, 107, 3902-3909. 35 VAMP 10.0, Dassault Systemes Biovia Corp. 2017. 36

Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB+, a Sparse Matrix-Based Implementation of

the DFTB Method J. Phys. Chem. A, 2007, 111, 5678-5684. 37 Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the Self-Consistent-Charge DensityFunctional Tight-Binding Method (SCC-DFTB) J. Chem. Theory Comput., 2011, 7, 931-948. 38 Brandenburg, J. G.; Hochheim, M.; Bredow, T.; Grimme, S. Low-Cost Quantum Chemical Methods for Noncovalent Interactions J. Phys. Chem. Lett., 2014, 5, 4275-4284. 39 Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB+, a Sparse Matrix-based Implementation of the DFTB Method. J. Phys. Chem. A, 2007, 111, 5678. 40

Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Electronic Structure Calculations on

Workstation Computers: The Program System Turbomole Chem. Phys. Lett., 1989, 162, 165169.


46

Page 47 of 49 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


41

TURBOMOLE

7.1

2017,

a

development

of

University

of

Karlsruhe

and

Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available form http://www.turbomole.com; COSMOlogic GmbH & Co. KG, Germany, 2017, accessed Dec 12 2017. 42 ANCOPT: Approximate Normal Coordinate Rational Function Optimization Program, S.Grimme, University of Bonn, 2005-13. 43 Eckert, F.; Pulay, P.; Werner, H.-J. Ab initio Geometry Optimization for Large Molecules J.Comput.Chem. 1997, 18, 1473-1483. 44 Not published; a brief description is given in the publication of GFN-xTB ref. 30. 45 Rauhut, G.; Clark, T.; Steinke, T. A Numerical Self-consistent Reaction Field (SCRF) Model for Ground and Excited States in NDDO-based Methods J. Am. Chem. Soc., 1993, 115, 9174-9181. 46 Klamt, A. Conductor-like Screening Model for Real Solvents: A New Approach to the Quantitative Calculation of Solvation Phenomena J. Phys. Chem., 1995, 99, 2224-2235. 47 Eckert, F.; Klamt, A. Fast Solvent Screening via Quantum Chemistry: COSMO-RS approach AIChE J., 2002, 48, 369-385. 48 Klamt, A. The COSMO and COSMO-RS Solvation Models WIREs: Comput. Mol. Sci., 2011, 1, 699-709. 49 COSMOlogic GmbH & Co. KG, Germany, 2017.


47


Page 48 of 49

50 Becke AD. Density-functional Exchange-energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A., 1988, 38, 3098. 51 Fraczkiewicz, R.; Lobell, M.; Göller, A. H.; Krenz, U.; Schoenneis, R.; Clark, R. D.; Hillisch, A. Best of Both Worlds: Combining Pharma Data and State of the Art Modeling Technology To Improve in Silico pKa Prediction J. Chem. Inf. Model., 2015, 55, 389-397. 52

53

ADMET Predictor, version 7.1; Simulations Plus, Inc.: Lancaster, CA, 2014. R version 3.4.0 (2017-04-21), Copyright (C) 2017 The R Foundation for Statistical

Computing, Platform: x86_64-redhat-linux-gnu (64-bit). 54 Goerigk, L.; Grimme, S. A Thorough Benchmark of Density Functional Methods for General Main Group Thermochemistry, Kinetics, and Noncovalent Interactions Phys. Chem. Chem. Phys., 2011, 13, 6670-6688. 55 Kesharwani M. K.; Karton, A.; Martin, J. M. L. Benchmark ab Initio Conformational Energies for the Proteinogenic Amino Acids through Explicitly Correlated Methods. Assessment of Density Functional Methods, J. Chem. Theory Comput., 2016, 12, 444–454. 56 Tao, J.; Perdew, J. P.; Staroverov, V. N.; Suseria, G. E. Phys. Rev. Lett., 2003, 91, 146401. 57 Labute, P. J. On the Perception of Molecules from 3D Atomic Coordinates J. Chem. Inf. Model., 2005, 45, 215-221. 58 Harder, E.; Damm, W.; Maple, J.; Wu, C.; Reboul, M.; Xiang, J. Y.; Wang, L.; Lupyan, D.; Dahlgren, M. K.; Knight, J. L.; Kaus, J. W.; Cerutti, D. S.; Krilov, G.; Jorgensen, W. L.; Abel,


48

Page 49 of 49 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


R.; Friesner, R. A. OPLS3: A Force Field Providing Broad Coverage of Drug-like Small Molecules and Proteins J. Chem. Theory Comput., 2016, 12, 281-296. 59 Storer, J. W.; Giesen, D. J.; Cramer, C. J.; Truhlar, D. G. Class IV Charge Models: A New Semiempirical Approach in Quantum Chemistry J. Comput.-Aided Mol. Des., 2004, 9, 87– 110. 60 Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-quality Atomic Charges. Am1-bcc Model: II. Parameterization and Validation J. Comput. Chem., 2002, 23, 1623– 1641. 61 Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces, B. Pullman, Ed., Reidel, Dordrecht, Holland, 1981, pp. 331-342.


49

Reliable and Performant Identification of Low-energy Conformers in

Recommend Documents