Evaluating Competing Intermolecular Interactions through Molecular

Article pubs.acs.org/crystal

Cite This: Cryst. Growth Des. XXXX, XXX, XXX−XXX

Evaluating Competing Intermolecular Interactions through Molecular Electrostatic Potentials and Hydrogen-Bond Propensities Published as part of a Crystal Growth and Design virtual special issue Honoring Prof. William Jones and His Contributions to Organic Solid-State Chemistry Bhupinder Sandhu,† Ann McLean,† Abhijeet S. Sinha,† John Desper,† Amy A. Sarjeant,‡ Shyam Vyas,‡ Susan M. Reutzel-Edens,§ and Christer B. Aakeröy*,† †

Department of Chemistry, Kansas State University, Manhattan, Kansas 66506, United States The Cambridge Crystallographic Data Centre, Piscataway, New Jersey 08854, United States § Small Molecule Design & Development, Eli Lilly and Company, Indianapolis, Indiana 46285, United States ‡

S Supporting Information *

ABSTRACT: The structural chemistry of six thiazole amides has been explored using experimental crystallographic data, and a combination of calculated hydrogen-bond propensities (HBPs) and hydrogen-bond energies. Both methods correctly identify the main hydrogen-bonded synthon, a pairwise N−H···N dimer, which appears in all the available structures. The strength and stability of the homosynthon in these compounds were tested by attempting to co-crystallize all six compounds with 20 different carboxylic acids. A total of 120 attempted reactions were carried out via liquid-assisted grinding experiments; sixty of these reactions produced a co-crystal, as indicated by IR spectroscopy. HBP calculations were employed in order to try to predict which of the 120 reactions would be successful. A multi-component score (MC score) was obtained from the hydrogen-bond propensity calculations, and this score coupled with a cut-off value >0.0 resulted in a 77% agreement between prediction and experiment (88% success for aliphatic acids and 67% for aromatic acids). By adjusting the MCcut‑off ≥ −0.1, the overall success rate for co-crystal prediction increased to 91% while significantly reducing the number of false negatives. The crystal structures of eight co-crystals were obtained, and all of them contain the synthon that was predicted by both hydrogen-bond energy and HBP calculations. The results demonstrate the importance of using complementary energy-based and structural chemistry methods for identifying which supramolecular pathways are accessible (and which are most likely) as small molecules transition from solution to the crystalline solid state.

1. INTRODUCTION It is well known that covalent synthetic modifications can be used to systematically modulate the electronic environment of a molecule, but when it comes to controlling and predicting the supramolecular assembly of molecules, the focus is on sitespecific reactivity and recognition of an intermolecular nature.1 In this context, supramolecular synthons2,3 have provided a foundation for the assembly of solid-state architectures with desired structural features, topologies, and chemical compositions.4−6 The hydrogen bond7−10 is, due to its strength and directionality, undoubtedly the most frequently used synthetic vector in supramolecular synthesis.11 Consequently, much work has been devoted to exploring how hydrogen-bond specificity is affected by geometric, steric, and electronic features within a particular molecular structure.12,13 However, it is still very challenging to predict which set of intermolecular interactions a given molecule is most likely to form in the solid state; how, in competition with other intermolecular forces, can hydrogen bonding be used to reliably direct molecular self-assembly in three dimensions? Without this knowledge, it will be extremely difficult to successfully tackle more complex questions such as © XXXX American Chemical Society

whether a molecule will crystallize or not, or if it is prone to polymorphism14,15 or not, issues which are of considerable importance to the manufacturing, delivery, and performance of any high-value molecular solids.16 Various approaches have been proposed to elucidate the structural outcome in cases where multiple intermolecular interactions are possible. Hunter et al. used a combination of experimentally derived parameters and electrostatic potential surfaces to predict supramolecular interactions in solid as well as in solution phase.17−19 Fábián employed molecular shape and polarity as benchmarks for predicting the formation of cocrystals.20,21 Price and co-workers focused on lattice energy comparisons of co-crystals and pure components as a basis for predicting co-crystal formation,22 and Velaga suggested that drugs and co-formers with similar Hanson solubility parameters are likely to form co-crystals with each other.23 Galek et al. demonstrated that hydrogen-bond propensity calculations,24,25 which rely on a statistical analysis of the occurrence of Received: October 17, 2017

A

DOI: 10.1021/acs.cgd.7b01458 Cryst. Growth Des. XXXX, XXX, XXX−XXX

Crystal Growth & Design

Article

Scheme 1. Representations of Three Postulated Synthons in a Generic Thiazole-Amidea

a (a) Synthon A, N−H (amide)···N(aromatic); (b) Synthon B, N−H (amide)···S(aromatic); and (c) Synthon C, N−H (amide)···OC(amide); (R1 = H, CH3; R2 = methyl, ethyl, or phenyl).

Scheme 2. Three Postulated Heteromeric Hydrogen-Bond Based Synthons between a Thiazole-Amide and a Carboxylic Acida

a

R1 = H, CH3; R2 = methyl, ethyl, or phenyl.

methods for predicting the likelihood of co-crystal formation and for rationalizing the observed intermolecular interactions. A structural informatics based approach to hydrogen-bond propensity calculations has been combined with a pairwise energy-based approach which relies on an electrostatic description of key hydrogen-bond interactions. The overall goal is to develop robust and transferable protocols for identifying and predicting which hydrogen bonds are most likely to appear in homomeric and heteromeric molecular solids when there are numerous potential avenues for assembly. First, we want to establish what intermolecular preference, if any, the N−H group displays and whether we can rationalize this in a convincing manner. The interaction between the hydrogen-bond donor and the three possible acceptors would lead to three different synthons, Scheme 1, taking into account only the trans CO-NH isomer and disregarding any catemeric versions of synthons A and B. We used structural informatics as a guideline to determine the most plausible conformation in these target molecules using data from the CSD. The search was restricted to structures without disorder, and with an Rfactor below 5%. According to our CSD search (CSD database 5.38, Nov 2016 with updates from Feb and May 2017),20,42−44 we found 6335 relevant structures containing the amide functionality, and only about 0.5% of these displayed a cis CO-NH isomer. Therefore, we felt justified in focusing exclusively on the trans CO-NH geometry in these target molecules. In addition, structural data complemented by geometry optimizations (see Section 2.2) indicated that the amide is likely to be co-planar with the thiazole. Second, we want to try to break the homomeric interactions (synthons A−C) by introducing a suitable carboxylic acid that

hydrogen bonds in relevant structures present in the Cambridge Structural Database (CSD),20,26−29 could be used for predicting the outcome of co-crystallizations of the drug lamotrigine,30,31 and Jones and co-workers employed hydrogenbond propensities as a way of quantifying homomeric and heteromeric interactions.32−34 Tilbury et al. recently used theoretical and data-driven propensity approaches to predict hydrates.35 Validation studies performed so far indicate that a combination of methods is valuable when conducting a knowledge based co-former screen.29 Thiazole rings are present in many pharmaceutical compounds, and drugs such as FB,36 AMG 517,37 meloxicam,38,39 and nitazoxanide,40,41 contain amide-substituted thiazoles as part of their chemical make-up. The multiple donor and acceptor atoms on the thiazole amide means that hydrogen-bond driven co-crystallizations can, in principle, follow different paths depending on which sites participate in the structure directing assemblies. In turn, this also means that multiple solid forms, with different physical properties, can be targeted in co-crystal synthesis. Knowing in advance which supramolecular pathways are accessible, and which of those possible options is the most likely route for a molecule, is thus clearly of importance for developing solid forms with optimum bulk properties. In the work presented herein we describe a systematic structural analysis of six multi-functional custom-designed thiazoles with three hydrogen-bond acceptors (amide carbonyl O, thiazole S, and thiazole N) and one hydrogen-bond donor (an amide N−H). The structural landscape of all these compounds is explored, and their ability to form co-crystals with a wide range of carboxylic acids has been evaluated. In parallel with the experimental work, we have employed two B



Article

Scheme 3. Six Thiazole Amides, T1−T6, (Top) and 20 Carboxylic Acid Co-Formers

targets (thiazole amides) and probes (carboxylic acids) is listed in Scheme 3. Molecules T1−T6 have systematically increasing molecular weight, and they also display specific and controllable changes to the charge density on the synthon. By exploring the structural landscape of this family of closely related molecules, we may be able to highlight the relative influence of small structural differences. Likewise, for the co-formers the aliphatic

may favor heteromeric interactions that can drive co-crystal formation (synthons D−F, Scheme 2). Our working strategy for the synthesis of co-crystals is based on optimizing the affinity between different molecules bearing complementary hydrogen-bonding moieties, thereby controlling the balance between homomeric (re-crystallization) and heteromeric (co-crystallization) outcomes. The complete set of C



Article

Table 1. Functional Groups Used To Determine the Hydrogen-Bond Propensities for the Six Target Moleculesa

The labels in the figures can be explained as follows: Tn = atom makes n bonds, c = atom is cyclic, @ = bond is acyclic, and Hn = n bonded hydrogen atoms. a

edited, functional groups were selected as suggested by Mercury (Table 1 represents functional groups for T1−T6, see Supporting Information for co-crystal functional groups), a training dataset (350600 structures per functional group) was made, and the propensity values were calculated using a logistic regression model with an ROC curve higher than 0.800 (“good discrimination”). The different amide substituents in targets T1/T2, T3/T4, and T5/T6 were modeled separately as shown in Table 1. The hydrogen-bond propensities for T1−T6 when screened against 20 different co-formers were calculated, and the multi-component score (MC score) was obtained by subtracting the propensity of the best homo-interaction (pure components, T:T or C:C) from the propensity of the best hetero-interaction (T:C or C:T), eq 4. A positive MC score indicates that the hetero-interaction needed for cocrystal formation is favored, whereas a negative MC score indicates that homo-interactions are favored and therefore hydrogen-bond directed co-crystallization is unlikely to occur.

carboxylic acids show increasing chain length and molecular flexibility. Similarly, the aromatic acids explore the addition of competing hydrogen bonding groups and the influence of substituent groups on the aromatic ring.

2. EXPERIMENTAL SECTION 2.1. Materials. 2-Amino-thiazole, 2-amino-5-methyl-thiazole, acetic anhydride, propionic anhydride, and benzoyl chloride were purchased from Aldrich and utilized without further purification. Synthetic procedures and characterization of all target molecules are provided in the Supporting Information. Dicarboxylic acids and benzoic acids (BA) were purchased from commercial sources and used as received. Melting points were measured using a Fisher-Johns melting point apparatus. Solution 1H NMR data were collected in DMSO-d6 on a Varian Unity plus 400 MHz NMR spectrometer. IR spectra of cocrystal screening experiments were recorded with a Nicolet 380 FT-IR spectrometer using an attenuated total reflection (ATR) technique and ZnSe as the crystal. 2.2. Molecular Electrostatic Potential Calculations. Molecular electrostatic potential surfaces (MEPs) of T1−T6 were generated with density functional B3LYP level of theory using 6-311++G** basis set in a vacuum. The initial geometry for each target molecule was selected after exploring the available energy space as a function of the rotatable exocyclic thiazole−nitrogen torsion angle. All molecules were then geometry optimized with the maxima and minima in the electrostatic potential surface (0.002 e/au isosurface) determined using a positive point charge in vacuum as a probe. The calculations yield the interaction energy (kJ/mol) between the positive point probe and the surface of the molecule at the point of contact. All calculations were carried out using Spartan’08 software (Wavefunction Inc.). 2.3. Use of Calculated Hydrogen-Bond Energies for Predicting Structural Outcomes. In order to determine which of the three postulated synthons is most likely to appear, based on interaction energies in the crystal structures of the pure compounds T1−T6, we calculated the α (hydrogen-bond donor) and β (hydrogen-bond acceptor) parameters from the maxima and minima on the MEPs, respectively (eqs 1 and 2), and the free energy of interaction is given by the product, −α β (eq 3).19 α = 0.0000162MEPmax 2 + 0.00962MEPmax

(1)

β = 0.000146MEPmin 2 − 0.00930MEPmin

(2)

E = − ∑ αiβj ij

⎛ propensity of the best ⎞ MC score = ⎜ ⎟ ⎝ heteromeric interaction ⎠ ⎛ propensity of the best ⎞ −⎜ ⎟ ⎝ homomeric interaction ⎠

(4)

2.5. Co-Crystal Screening and X-ray Crystallography. T1−T6 were put through a co-crystal screen using liquid-assisted grinding (LAG) in a drop of methanol against 10 aliphatic dicarboxylic acids and 10 aromatic acids (Scheme 3). In all 120 reactions the reactants were mixed in stoichiometric ratios (2:1 for diacids and 1:1 for monoacids), and the solid resulting from each reaction was characterized using IR spectroscopy to determine if a co-crystal had formed. The IR analysis focused on the positions of readily identifiable vibrational modes (cf. OC (amide), OC (acid), etc.) of the pure thiazole amides and acids relative to the corresponding bands of the resulting solids. Band shifts greater than seven wavenumbers as well as the appearance of broad stretches around 2300 and 1800 cm−1 as a result of intermolecular O−H···N hydrogen bonds were taken as an indication of co-crystal formation (such hydrogen bonds would not be possible in either of the pure compounds, see Supporting Information). Subsequently, all 120 combinations from the grinding experiments were dissolved in a minimum amount of methanol or nitromethane (1−2 ml, both reactants were readily soluble in these solvents) and kept in a vial for slow evaporation in order to obtain crystals suitable for single crystal X-ray diffraction. Of the 120 LAG reactions, 60 produced co-crystals according to IR data and eight of those, yielded crystals suitable for single-crystal X-ray diffraction. X-ray experimental data and crystallographic data are given in the Supporting Information.

(3)

The same approach was used to evaluate the energies of heteromeric synthons. The interaction energies in the co-crystal were determined by pairing the best donor of carboxylic acid with the best acceptor of the target molecule and best acceptor of carboxylic acid with the best donor of the target molecule. 2.4. Use of Hydrogen-Bond Propensities for Predicting Structural Outcomes. To complement the electrostatic-based calculations, we also employed hydrogen-bond propensity24,29,30 calculations (CSD Version 5.38 and Mercury 3.9) as a way of predicting the most likely interactions in the structures of the parent molecules, T1−T6, and their co-crystals with the 20 hydrogen-bond donors shown in Scheme 3. Each compound was sketched and auto-

3. RESULTS AND DISCUSSION First, we wanted to elucidate if the six target molecules themselves favor a particular synthon in their solid-state structures, and if so, could those interactions be predicted or rationalized using an electrostatic or propensity-based approach? Second, we wanted to establish if we can break the homomeric interactions between target molecules and form cocrystals with aliphatic or aromatic carboxylic acids, and if so, D



Article

synthon A is more likely to occur in these target molecules due to the simultaneous involvement of two hydrogen bonds (also referred to as the “supramolecular chelating effect”37). 3.1.3. Hydrogen-Bond Propensities. The propensity calculations for the six target molecules consider all possible pairwise interactions that may take place between one of the three acceptors and the sole hydrogen-bond donor, N−H (amide). The results show that the individual hydrogen bonds that make up synthon A have the highest propensity followed by those in synthon C and then in synthon B, Table 2. Therefore, the structural informatics approach is consistent with the energy calculations in suggesting that self-complementary synthon A is favored over the postulated alternatives in the crystal structures of T1−T6. The differences in propensity values of a given synthon across T1−T6 show the sometimes large effect that pendant groups can have. The knowledge-based approach that serves as the foundation for the propensity values also includes parameters such as “steric density”, that would not be easy to identify and evaluate without an automated and systematic analysis of thousands of structural data points. 3.1.4. Experimentally Observed Crystal Structures. We obtained crystal structures for T2, T4, and T6 (crystal structures for T1 (YODJAD)45 and T5 (NORLAI)46 have previously been reported) and in all five cases, synthon A (with a graph-set notation of R22(8)) was observed, Figure 2. The only conventional hydrogen-bond donor in this group of molecules, the N−H amide moiety, selects the N(thiazole) acceptor site as its preferred partner. The preference for synthon A over synthon B is readily explained on the basis of electrostatics. As there is no geometric bias for either motif, the outcome appears to be determined by the fact that the sulfur acceptor atom with its significantly lower (less negative) electrostatic potential value is not competitive relative to the thiazole ring nitrogen atom. This is also borne out in the propensity analysis which shows a very low probability for the occurrence of the hydrogen bonds in synthon B, implying that it is not as frequently observed in crystal structures. Even though individual N−H···N interactions (in synthon A) are weaker than individual N−H···OC hydrogen bonds (in synthon C), synthon A gains an advantage as a result of the geometry and relative position of the donors/ acceptors on these molecules which facilitates the appearance of two hydrogen bonds through a supramolecular “chelating” effect. We hypothesize that such dimeric assemblies can most likely compensate for individually weaker interactions in a substantial number of cases, but only if the electrostatic differences between the possible acceptors is not too large. It is also possible that this energetic bias is most important during crystal growth, where a dimeric synthon A growth unit would be preferred over synthon C. It is, of course, not too surprising that a pair of weak interactions can overcome a single interaction that is thermodynamically more favorable. However, it is a real challenge to predict exactly where the energetic “tipping point” may reside, and this is why a combined MEPs/ HBP approach is potentially very powerful. If/when both methods agree (i.e., they rank the order of interactions in the same manner), this likely means that combinations of interactions have limited ability to overcome the stabilization of isolated stronger interactions. If the HBP shows that weaker interactions are represented in the statistics, it would indicate that crystal packing forces can tilt the balance in favor of the weaker interactions.

could those new synthons be predicted or rationalized using an energy or propensity-based approach? The following sections will deal with each protocol separately as they were applied to co-crystals of T1−T6 with aliphatic diacids and aromatic monoacids, respectively. 3.1. Structural Chemistry and Interaction Analyses of T1−T6. 3.1.1. Molecular Electrostatic Potentials. In order to obtain the necessary molecular electrostatic surface values for T1−T6, we performed energy minimizations on the most likely conformations (as determined from an analysis of relevant data from the CSD) of these target molecules with the restriction that all have a trans-configuration of the amide group (as suggested by the CSD analysis). The most stable conformation for each molecule was subsequently used for the MEPs calculations (relevant energies tabulated in the Supporting Information). It is worth noting that the energy optimized conformations are not necessarily completely identical to those that may appear in the solid state, where a variety of close contacts and packing forces may distort some geometric parameters away from idealized gas phase values. However, these idealized conformations are likely to be most relevant in the solution phase at the point when target molecules begin to recognize and bind to each other during nucleation and growth. MEPs values for each acceptor group (aromatic N, aromatic S and OC) and the one N−H donor are provided in Figure 1. According to Etter’s guidelines, the best acceptor will bind to

Figure 1. Electrostatic potentials values (in kJ/mol); (a) T1, (b) T2, (c) T3, (d) T4, (e) T5, and (f) T6.

the best donor in a competitive environment, and if we adopt an electrostatically-based ranking of the three acceptors, the outcome is as follows: CO > N(aromatic) ≫ S(aromatic). Therefore, in a system without additional geometric bias, the most likely hydrogen bond would be the N−H (amide)···O C interaction (synthon C, Scheme 1). 3.1.2. Hydrogen-Bond Energies. The calculated hydrogenbond energy for individual synthons in T1−T6 is tabulated in the Supporting Information. The average hydrogen-bond energy of synthon A across T1−T6 is −29.75± 1.17 kJ/mol (−14.88 ± 0.58 kJ/mol for each N−H···N bond), synthon B is −10.76 ± 0.67 kJ/mol (−5.38 ± 0.34 kJ/mol for each N−H···S bond), and synthon C is −19.41 ± 1.47 kJ/mol, Table 2. Although the hydrogen-bond energy of N−H···OC is greater than each of the individual bonds within synthons A and B, E



Article

Table 2. Hydrogen-Bond Energies (Higher Absolute Value Indicates Stronger Interaction) and Propensities (Larger Value Indicates Increased Likelihood of Formation) for Each Postulated Synthon in T1−T6a

a

All propensities are given for individual hydrogen bonds whereas hydrogen-bond energies are given per synthon.

Figure 2. Synthon A appears in the crystal structures of T1−T2, T4−T6.

Scheme 4. Likely (D and F) and Unlikely (E) Synthons in Co-Crystals of T1−T6

3.2. Structural Chemistry and Interaction Analyses of Co-Crystals of T1−T6. Having shown the statistically and energetically favored synthon A exclusively formed in T1/T2/ T4/T5/T6, we wanted to determine if this homomeric interaction could be disrupted by a sufficiently competitive co-former to drive the synthesis of co-crystals of T1−T6 under the same experimental conditions (i.e., liquid assisted grinding method with methanol). We explored this idea through a

systematic co-crystal screen against 20 carboxylic acids, initially postulating that three different motifs could result from cocrystallizations between a carboxylic acid and T1−T6, but given the behavior of the target molecules themselves, we were able to rule out the likelihood of finding synthon E in a co-crystal, Scheme 4. A heterodimer, synthon D, could form by N− H(amide)···OC(acid) and O−H(acid)···N(thiazole) hydrogen bonds. Alternatively, O−H(acid)···OC(amide) interF



Article

Table 3. Hydrogen-Bond Energies of Each Hetero-Synthon in T1−T6 Co-Crystals with Aliphatic and Aromatic Acids co-former aliphatic acids aromatic acids hydroxybenzoic acids aminobenzoic acids

Synthon D (kJ/mol) −36.94 −33.04 −33.85 −35.03

± ± ± ±

0.73 0.67 0.67 0.72

Synthon F (kJ/mol) −24.51 −25.06 −22.26 −20.35

± ± ± ±

1.10 1.11 0.99 0.91

Synthon G/I (kJ/mol)

Synthon H/J (kJ/mol)

N/A N/A −25.35 ± 1.13 −17.23 ± 0.76

N/A N/A −19.47 ± 0.90 −13.23 ± 0.61

Scheme 5. Possible Synthon G and H (in hydroxyBA) and Synthon I and J (in aminoBA) in Co-Crystals of T1−T6

Table 4. Attempted Co-Crystallizations Using LAG (Methanol) of T1-T6 with Aliphatic and Aromatic Acids

G



Article

Figure 3. Primary hydrogen bonds in the crystal structures of (a) T6-Suc, (b) T6-Sub, (c) T6-Seb, (d) T6-Aze, and (e) T6-Dod.

produce co-crystals with aliphatic acids under the standard conditions employed here. The experimental screening of the 10 aromatic acids against T1−T6 yielded 38 co-crystals from 60 experiments, Table 4. The success rate for the aromatic acids (63%) is higher than that for the aliphatic acids (37%) probably due to the presence of additional hydrogen-bond donors in four of them which generally increases the chance of hydrogen bonding. Moreover, all target molecules T1−T6 show better success with aromatic acids, particularly with stronger aromatic acids with lower pKa values (and larger positive electrostatic potential values on their acidic protons). Overall, 60 out of 120 experiments yielded positive outcomes, a rather high success rate considering the very limited scope of the experimental screen, only LAG from methanol was used to screen for co-crystals and only one solvent (methanol or nitromethane) was used to grow single crystals. With further experimentation, it is reasonable to assume that some pairs that did not succeed in co-crystallizing in this screen could eventually form, and that polymorphs of any of the pairs might be found. Still, it is necessary to employ consistent experimental conditions when trying to elucidate the success or suitability of a protocol for co-crystal prediction. 3.2.3. Experimentally Observed Crystal Structures: CoCrystals with Aliphatic Diacids. We obtained single-crystal XRD data for five co-crystals with aliphatic acids, and in all five crystal structures the carboxylic acid O−H···N thiazole and amide N−H···OC acid moieties give rise to the robust R22(8) hydrogen-bonded motif in a 2:1 target/acid stoichiometry, Figure 3. The carboxylic groups of each co-former provide ideal hydrogen-bond complementarity to the thiazole-amide, and the combination of O−H···N and N−H···O hydrogen bonds in the heterodimer favor synthon D over synthon F. 3.2.4. Experimentally Observed Crystal Structures: CoCrystals with Aromatic Acids. Three crystal structures were obtained with aromatic acids with varying stoichiometry (1:2 ratio for T1-3HydroxyBA, 2:1 for T2-4HydroxyBA, and 1:1 for T4-3NitroBA). As with the aliphatic acid co-crystals, a robust R22(8) hydrogen-bonded synthon D is formed using O−H···N

actions could form synthon F (while keeping the synthon A homodimer intact). 3.2.1. Hydrogen-Bond Energies for Co-Crystals. The hydrogen-bond energy for synthon D and synthon F in each co-crystal combination (120 combinations) was calculated and is tabulated in the Supporting Information (only heteromeric interactions are considered). The more negative the hydrogenbond energy, the stronger the interaction between the two different components, favoring co-crystallization. For the aliphatic acids, the energy of synthon D is −36.94 ± 0.73 kJ/ mol (−16.82 ± 0.84 kJ/mol for N−H···OC bond and −17.87 ± 0.83 kJ/mol for O−H···N bond), and for heterosynthon F it is −24.51 ± 1.10 kJ/mol, Table 3. For the aromatic acids without an additional donor (−OH or −NH2 substituent) the energy of synthon D is −33.04 ± 0.67 kJ/mol and for synthon F is −25.06 ± 1.11 kJ/mol. Although the individual hydrogen bonds comprising synthon D are weaker than the O−H(acid)···OC(amide) bond in synthon F, the combination of these hydrogen bonds in synthon D is more stabilizing than the single interaction in synthon F.37 In order to explore how introducing other competitive hydrogen-bonding donor groups could tilt the balance of interactions favoring (or disfavoring) co-crystal formation, we also examined T1−T6 in combination with hydroxy and aminobenzoic acids. Competitive hydrogen-bonding arrangements, synthons G and H for hydroxyBA and synthons I and J for aminoBA, are shown in Scheme 5. The hydrogen-bond energies were calculated for D, F, G, and H (for hydroxyBA) and D, F, I, and J (for aminoBA) and are tabulated in the Supporting Information. As summarized in Table 3, synthon D is still the strongest interaction suggesting that the presence of hydroxyl and amino groups is not likely to alter the ranking. 3.2.2. Experimental Co-Crystal Screening of T1−T6 against 10 Aliphatic and 10 Aromatic Carboxylic Acids. The experimental screen produced co-crystals in 22 of the 60 experiments between T1−T6 and aliphatic diacids, Table 4 (from IR). T5 and T6 readily co-crystallized with 90% and 100 % success rates, respectively. By contrast, T3/T4 did not H



Article

Figure 4. Main hydrogen bonds in the crystal structures of (a) T1-3HydroxyBA, (b) T2-4HydroxyBA, and (c) T4-3NitroBA.

Table 5. Summary of Predicted vs Experimental Synthons in Co-Crystals of T1−T6 with Five Aliphatic and Three Aromatic Acids T1−T6

hydrogen-bond energies

hydrogen-bond propensities

experimental

no. of experimental data points

aliphatic acids aromatic acids hydroxybenzoic acids

Synthon D Synthon D Synthon D

Synthon D Synthon D Synthon D

Synthon D Synthon D Synthon D and Synthon G

5 1 2

Figure 5. Confusion matrices determined from multi-component HBP results for T1−T6 molecules with co-formers, (a) explanation of matrix (b) MCcutoff: > 0.00, (c) MCcutoff: ≥ −0.10.

the crystal structure of T2-4HydroxyBA, both synthons D and synthon G were observed. 3.2.5. Comparison of Synthon Energies of Co-Crystals with Aliphatic Acids vs. Aromatic Acids. For aliphatic acids, the dimeric synthon D was favored over the monomeric synthon F by 13 kJ/mol. In five crystal structures obtained with T6, synthon D was observed irrespective of the aliphatic acid.

and N−H···O hydrogen bonds in each case (Figure 4). In T13HydroxyBA, the additional hydrogen-bond donor on 3HydroxyBA binds to the carbonyl oxygen of a second acid, OH(hydroxyl)···OC(acid), Figure 4a. In T2-4HydroxyBA, the extra hydrogen-bond donor interacts with a carbonyl oxygen atom of a second T2 to form 1-D chains, Figure 4b. In I



Article

Figure 6. Multi-component HBP analysis results for a co-crystal screen of T1−T6 against aliphatic acids; true positive and true negatives (green triangles); false negatives (yellow triangles); and false positives (red triangles). More details on the MC-HBP results are given in the Supporting Information.

prediction but positive via experimental), and true negative (TN; negative for co-crystal via both prediction and experimental), Figure 5a. The MC score was divided into two categories; MCcutoff > 0.00 and MCcutoff ≥ −0.10. Any MC value >0.00 was considered as a prediction of co-crystal formation and an MC value below or equal to 0.00 was considered as a prediction against co-crystal formation. In 93 out of the 120 combinations (77%), the HBP model accurately predicted the outcome of the reaction; 41 out of 93 combinations were true positive and 52 of the 93 combinations were true negative, respectively, Figure 5b. Experimentally, T1−T4 displayed a 7.5% and T5−T6 a 95% success rate of co-crystal formation with aliphatic acids. Such a discrepancy between molecules that carry the same principle functionality is difficult to predict or foresee using chemical intuition, but we were delighted that the HBP model delivered predictions that reflected these highly unexpected experimental results very well. Using an MCcutoff > 0.00, only 1 out of 40 reactions (2.5%) between an aliphatic diacid and T1−T4 was predicted to yield a co-crystal, whereas 18 of 20 combinations (90%) between a diacid and T5−T6 were predicted to produce a co-crystal. There was an 88% agreement (53 out of 60 combinations were either true positive or true negative) between the predicted and experimental results for aliphatic acids. Among the seven unsuccessful predictions, two were false positives (red triangles) and five were false negatives (yellow triangles), Figure 6. In T1−T4 with aliphatic acids, homomeric interactions; N−H(amide)···N(thiazole) or O−H(acid)···O C(acid) are more probable than the heteromeric interaction, O−H(acid)···N(thiazole), whereas in T5−T6, the heteromeric interaction, O−H(acid)···N(thiazole), dominates over homomeric interactions which suggests that descriptors such as “steric density” and “aromaticity” play a significant role in the co-crystallization outcome. The influence of these features on the different assembly paths that are available for potential cocrystals of T1−T6 is extremely difficult to quantify and evaluate using computational methods. However, this knowledge-based HBP protocol has allowed us to accurately bring out the experimental consequences of very subtle changes in molecular structure. The experimental co-crystal screening of 10 aromatic acids against T1−T6 produced 38 positive results out of 60 attempts.

For aromatic acids with only one hydrogen-bond donor and acceptor, synthon D was still favored over synthon F by 8 kJ/ mol and in the T3-3NitroBA co-crystal, only synthon D was observed. On the basis of hydrogen-bond energies, synthon D is favored for both hydroxybenzoic acid co-crystals (over synthons G, F, and H) and aminobenzoic acid co-crystals (over synthons F, I, and J). Two crystal structures of hydroxybenzoic acid co-crystals were obtained; in T1-3HydroxyBA, synthon D was observed as predicted. All conventional hydrogen-bond donors and acceptors are involved in heteromeric interactions (thus following Etter’s rule) in T2-4HydroxyBA. Without geometric constraints, the best donor of the acid (hydroxylic O−H) binds to the best acceptor of T2 (CO amide), the second-best donor of the acid (carboxylic O−H) binds to the second best acceptor of T2 (aromatic N), and the best donor of T2 (N−H amide) binds to the best acceptor of the acid (O C). Overall, the predictions made using the MEPs approach agree very well with the experimentally observed crystal structures, Table 5. 3.2.6. Evaluating Propensity Predictions for Co-Crystals with Aliphatic and Aromatic Acids. Hydrogen-bond propensity calculations were performed on T1-T6 against 20 coformers. The two highest propensity donor−acceptor interactions were taken into account for homomeric systems: N− H···N for the target molecule designated as T:T (target donor/ target acceptor) and O−H···OC or N−H···OC (in the case of aminoBA) for the acid co-former designated as C:C (coformer donor/co-former acceptor). The interaction with the maximum propensity for heteromeric systems was designated as T:C (target donor/co-former acceptor) or C:T (co-former donor/target acceptor). An MC score (as described in section 2.4) is used to predict whether a reaction (in this case, cocrystal formation) will or will not take place when two reactants are combined. A summary of the multi-component hydrogenbond propensity (HBP) screening results for T1−T6 with 20 co-formers is given in the Supporting Information. The predictions made using the HBP model were compared to experimental co-crystal screening results. A confusion matrix was used to analyze the data, as reported by Wood et al.29 The entries in the matrix are labelled true positive (TP; positive for co-crystal via both prediction and experimental), false positive (FP; positive for co-crystal via prediction but negative via experimental), false negative (FN; negative for co-crystal via J



Article

Figure 7. Multi-component HBP analysis results for a co-crystal screen of T1−T6 against aromatic acids; true positive and true negatives (green diamonds); false positives (red diamonds); and false negatives (yellow diamonds). More details on the MC-HBP results are given in the Supporting Information.

On the basis of hydrogen-bond propensity calculations, it was found that the heteromeric O−H(acid)···N(thiazole) hydrogen bond plays a more important role in deciding the outcome of an attempted co-crystallization than does the N− H(amide)···OC(acid) hydrogen bond. Although HBP values, which are based on probabilities, cannot be directly linked to thermodynamic parameters, it is reasonable, given the evidence presented herein, to assume that the former hydrogen bonds provide more of an enthalpic advantage than the latter and this may be translated into intermolecular bond strength, or binding efficiency. This observation may contribute to the development of more versatile guidelines for designing and selecting co-formers that are optimized for binding to a specific target molecule. The HBP results also show that synthon D is favored over synthon F in all aliphatic acid co-crystals and in co-crystals with aromatic acids that do not have additional donors, because the propensity of O−H(acid)···N (thiazole) is higher than that of the O−H(acid)···OC(amide) hydrogen bond in each cocrystal combination. The propensity results also favor D synthon over F/G/H, which are possible in hydroxyBA, and synthon J is favored over synthon D/F/I, which are possible in aminoBA co-crystals. The HBP protocol is in generally good agreement with the electrostatically based hydrogen-bond energy method; both agree that synthon D is favored over other possible synthons and this was supported by single crystal X-ray diffraction studies. The one instance where the two methods did not agree was in the case of T1−T6 and aminobenzoic acids; synthon D was predicted using the energy-based approach and synthon J was favored using the HBP model. Unfortunately, we have not been able to grow crystals suitable for single-crystal X-ray diffraction so we do not yet have the data to resolve this pending issue. Overall, the HBP model offers very encouraging results and when it is combined with electrostatic-based energy calculations on the most important homomeric and heteromeric interactions, the confidence in the predictions is significantly increased. Moreover, in those cases where crystallographic data were available for new co-crystals of T1−T6, the specific functional group interactions and synthons predicted by both methods were consistently observed. The consistency and correlation between hydrogen-bond energies

T5−T6 showed a higher success rate than T1−T4, but almost all target molecules were able to form co-crystals with a success rate greater than 50% (the one exception, T3, had a 40% success rate). The HBP models predicted mostly negative outcomes for T1−T4, whereas reactions of T5−T6 were predicted as positive for co-crystal formation (3AminoBA and 4AminoBA co-formers were always predicted to produce cocrystals, regardless of target). Overall, the outcome of 40 of the 60 reactions (67%) was accurately predicted by the HBP model (either true positives or true negatives) using an MC > 0.00 cutoff. Additionally, 16 of the 60 aromatic acid/T1−T6 combinations were false negatives (yellow diamonds) and four (red diamonds) were false positives (Figure 7). As indicated by Wood et al.,29 false negatives are a concern with prediction tools when it comes to creating a list of potential coformers for a systematic search for new solid forms, as these coformers will not be included in the experimental study. On the other hand, a false positive can also be detrimental as extensive experimental resources and time may then be wasted on a particular co-former that is expected to produce a co-crystal, but does not. Although the predictions suggest that the pairs are favored to form co-crystals, they do not indicate how much effort is required to find them. The choice of cut-off value for the MC score will clearly influence the prediction accuracy, that is the number of accurately predicted outcomes and the number of false positives and negatives. Using an MC score greater than 0.0, the HBP model gave 17% false negatives (21/120), and all but five of these reactions included an aromatic acid, Figure 5b. In order to avoid removing suitable co-formers that do form co-crystals but were excluded based on the initially selected 0.0 cut-off value, we noted that a modified value of MCcutoff ≥ −0.10 would reduce the number of false negatives from 21 to 5 and improve the overall success of the HBP predictions to 91% (109 out of 120 reactions), Figure 5c. However, we recognize that the MCcutoff value that offers the best prediction cannot be known a priori and is likely to be different for other pairs of molecules. This means that the value of the HBP approach lies in the relative ranking of co-formers, as opposed to the absolute MC scores. Such rankings could be used to prioritize experimental screens, especially in early-phase drug discovery when only small amounts of materials are available. K



Article

and the HBP protocol as guidance for predicting co-crystal formation indicate that the prevalence and efficiency of synthons that can drive the synthesis of multicomponent solids is strongly influenced by hydrogen-bond strength between cocrystal components.

Experimental data, IR data, hydrogen-bond energies data, hydrogen-bond propensities data (PDF) Accession Codes

CCDC 1578893−1578903 contain the supplementary crystallographic data for this paper. These data can be obtained free of charge via www.ccdc.cam.ac.uk/data_request/cif, or by emailing [email protected], or by contacting The Cambridge Crystallographic Data Centre, 12 Union Road, Cambridge CB2 1EZ, UK; fax: +44 1223 336033.

4. CONCLUSIONS The structural chemistry of six target molecules has been analyzed using hydrogen-bond energies and hydrogen-bond propensities as a way of predicting which intermolecular interactions are likely to appear in the solid state. On the basis of hydrogen-bond energies, synthon A was predicted to be favored over synthons B and C, thanks to a combination of enthalpic factors and a supramolecular chelating effect. The HBP results were consistent with the energy calculations, and offered a ranking of the hydrogen-bond donor/acceptor pairs suggesting that the synthons to which they belong can be ranked as follows: A > C > B. The two methods accurately predicted the experimentally observed intermolecular hydrogen-bond interactions in all five available crystal structures; only synthon A was present. We attempted to co-crystallize the six target molecules with 20 different carboxylic acids to probe the strength and stability of the homomeric interactions. The HBP calculations also predicted synthon D to be more favorable in all but the aminoBA co-crystals. Synthon J is favored over synthons D/F/I in aminoBA co-crystals. The HBP results were contrasted to the results of the experimental co-crystal screen comprising six target molecules and 20 carboxylic acids. An MC score produced from the HBP calculations using a cut-off value of >0.0 resulted in a 77% agreement between prediction and experiment (88% success for aliphatic acid and 67% for aromatic acids). By adjusting the MCcut‑off to ≥ −0.1, the overall success rate for co-crystal prediction increased to 91% while significantly reducing the number of false negatives. Eight crystal structures of T1−T6 co-crystals were obtained, each containing synthon D as the primary hydrogen-bond interactions; this structural consistency was predicted by both hydrogen-bond energy and HBP calculations. The structure determinations also reveal that each co-former provides an effective complement to the main molecular-recognition sites on T1−T6. The ability to make informed predictions about how different molecules recognize and bind to each other, and how they subsequently assemble into solid-state architectures is of critical importance in areas such as drug design and formulation, and the recently published guidance by the FDA on the regulatory classification of pharmaceutical co-crystals47 will continue to keep this fundamental science in sharp focus. Therefore, the successful use of structural informatics tools such as hydrogenbond propensity and calculated molecular electrostatic potential surfaces for mapping out the structural landscape of these types of molecules will have significant practical applications.

■

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Susan M. Reutzel-Edens: 0000-0003-0806-5565 Christer B. Aakeröy: 0000-0002-8947-2068 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We are grateful for financial support from Eli Lilly and Company and the Cambridge Crystallographic Data Centre (Contract number: 34397).

■

REFERENCES

(1) Aakeröy, C. B.; Epa, K. Springer: Berlin, Heidelberg, 2014; pp 125−147. (2) Shattock, T. R.; Arora, K. K.; Vishweshwar, P.; Zaworotko, M. J. Cryst. Growth Des. 2008, 8, 4533−4545. (3) Desiraju, G. R. Angew. Chem., Int. Ed. Engl. 1995, 34, 2311−2327. (4) Nangia, A.; Desiraju, G. R. Springer: Berlin, Heidelberg, 1998; pp 57−95. (5) Thalladi, V. R.; Goud, B. S.; Hoy, V. J.; Allen, F. H.; Howard, J. A. K.; Desiraju, G. R. Chem. Commun. 1996, 3, 401−402. (6) Reddy, D. S.; Craig, D. C.; Desiraju, G. R. J. Am. Chem. Soc. 1996, 118, 4090−4093. (7) Lehn, J. M. Science 2002, 295, 2400−2403. (8) Smith, G.; Baldry, K. E.; Byriel, K. A.; Kennard, C. H. L. Aust. J. Chem. 1997, 50, 727−736. (9) Bosch, E. CrystEngComm 2007, 9, 191−198. (10) Buckingham, A. D.; Del Bene, J. E.; McDowell, S. A. C. Chem. Phys. Lett. 2008, 463, 1−10. (11) Merz, K.; Vasylyeva, V. CrystEngComm 2010, 12, 3989−4002. (12) Bedeković, N.; Stilinović, V.; Piteša, T. Cryst. Growth Des. 2017, 17, 5732−5743. (13) Lehn, J.M. Rep. Prog. Phys. 2004, 67, 249−265. (14) Moulton, B.; Zaworotko, M. J. Chem. Rev. 2001, 101, 1629− 1658. (15) Aitipamula, S.; Nangia, A. In Supramolecular Chemistry; John Wiley & Sons, Ltd., 2012. (16) Nauha, E.; Bernstein, J. Cryst. Growth Des. 2014, 14, 4364− 4370. (17) Hunter, C. A. Angew. Chem., Int. Ed. 2004, 43, 5310−5324. (18) McKenzie, J.; Feeder, N.; Hunter, C. A. CrystEngComm 2016, 18, 394−397. (19) Musumeci, D.; Hunter, C. A.; Prohens, R.; Scuderi, S.; McCabe, J. F. Chem. Sci. 2011, 2, 883−890. (20) Macrae, C. F.; Bruno, I. J.; Chisholm, J. A.; Edgington, P. R.; McCabe, P.; Pidcock, E.; Rodriguez-Monge, L.; Taylor, R.; van de Streek, J.; Wood, P. A. J. Appl. Crystallogr. 2008, 41, 466−470. (21) Fábián, L. Cryst. Growth Des. 2009, 9, 1436−1443. (22) Issa, N.; Karamertzanis, P. G.; Welch, G. W. A.; Price, S. L. Cryst. Growth Des. 2009, 9, 442−453. (23) Mohammad, M. A.; Alhalaweh, A.; Velaga, S. P. Int. J. Pharm. 2011, 407, 63−71.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.cgd.7b01458. L



Article

(24) Galek, P. T. A.; Allen, F. H.; Fábián, L.; Feeder, N. CrystEngComm 2009, 11, 2634−2639. (25) Chemburkar, S. R.; Bauer, J.; Deming, K.; Spiwek, H.; Patel, K.; Morris, J.; Henry, R.; Spanton, S.; Dziki, W.; Porter, W.; Quick, J.; Bauer, P.; Donaubauer, J.; Narayanan, B. A.; Soldani, M.; Riley, D.; McFarland, K. Org. Process Res. Dev. 2000, 4, 413−417. (26) Allen, F. Acta Crystallogr., Sect. B: Struct. Sci. 2002, 58, 380−388. (27) Mapp, L. K.; Coles, S. J.; Aitipamula, S. Cryst. Growth Des. 2017, 17, 163−174. (28) Vologzhanina, A. V.; Sokolov, A. V.; Purygin, P. P.; Zolotarev, P. N.; Blatov, V. A. Cryst. Growth Des. 2016, 16, 6354−6362. (29) Wood, P. A.; Feeder, N.; Furlow, M.; Galek, P. T. A.; Groom, C. R.; Pidcock, E. CrystEngComm 2014, 16, 5839−5848. (30) Galek, P. T. A.; Chisholm, J. A.; Pidcock, E.; Wood, P. A. Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater. 2014, 70, 91−105. (31) Galek, P. T. A.; Pidcock, E.; Wood, P. A.; Bruno, I. J.; Groom, C. R. CrystEngComm 2012, 14, 2391−2403. (32) Delori, A.; Galek, P. T. A.; Pidcock, E.; Patni, M.; Jones, W. CrystEngComm 2013, 15, 2916−2928. (33) Delori, A.; Galek, P. T. A.; Pidcock, E.; Jones, W. Chem. - Eur. J. 2012, 18, 6835−6846. (34) Eddleston, M. D.; Arhangelskis, M.; Fábián, L.; Tizzard, G. J.; Coles, S. J.; Jones, W. Cryst. Growth Des. 2016, 16, 51−58. (35) Tilbury, C. J.; Chen, J.; Mattei, A.; Chen, S.; Sheikh, A. Y. Cryst. Growth Des. 2017, DOI: 10.1021/acs.cgd.7b00517. (36) Peterson, M. L.; Stanton, M. K.; Kelly, R. C.; Staples, R.; Cheng, A. CrystEngComm 2011, 13, 1170−1180. (37) Stanton, M. K.; Bak, A. Cryst. Growth Des. 2008, 8, 3856−3862. (38) Cheney, M. L.; Weyna, D. R.; Shan, N.; Hanna, M.; Wojtas, L.; Zaworotko, M. J. Cryst. Growth Des. 2010, 10, 4401−4413. (39) Tumanov, N. A.; Myz, S. A.; Shakhtshneider, T. P.; Boldyreva, E. V. CrystEngComm 2012, 14, 305−313. (40) Suresh, K.; Mannava, M. K. C.; Nangia, A. Chem. Commun. 2016, 52, 4223−4226. (41) Félix-Sonda, B. C.; Rivera-Islas, J.; Herrera-Ruiz, D.; MoralesRojas, H.; Höpfl, H. Cryst. Growth Des. 2014, 14, 1086−1102. (42) Taylor, R.; Macrae, C. F. Acta Crystallogr., Sect. B: Struct. Sci. 2001, 57, 815−827. (43) Bruno, I. J.; Cole, J. C.; Edgington, P. R.; Kessler, M.; Macrae, C. F.; McCabe, P.; Pearson, J.; Taylor, R. Acta Crystallogr., Sect. B: Struct. Sci. 2002, 58, 389−397. (44) Macrae, C. F.; Edgington, P. R.; McCabe, P.; Pidcock, E.; Shields, G. P.; Taylor, R.; Towler, M.; van de Streek, J. J. Appl. Crystallogr. 2006, 39, 453−457. (45) Yunus, U.; Tahir, M. K.; Bhatti, M. H.; Wong, W.Y. Acta Crystallogr., Sect. E: Struct. Rep. Online 2008, 64, o1516. (46) Zonouzi, A.; Mirzazadeh, R.; Rahmani, H.; Ng, S. W. Acta Crystallogr., Sect. E: Struct. Rep. Online 2009, 65, o817. (47) Gadade, D. D.; Pekamwar, S. S. Adv. Pharm. Bull. 2016, 6, 479− 494.

M


Evaluating Competing Intermolecular Interactions through Molecular

Recommend Documents