Rational Coformer or Solvent Selection for Pharmaceutical Cocrystallization or Desolvation

It is demonstrated that the fluid-phase thermodynamics theory conductor-like screening model for real solvents (COSMO-RS) as implemented in the COSMOtherm software can be used for accurate and efficient screening of coformers for active pharmaceutical ingredient (API) cocrystallization. The excess enthalpy, Hex, between an API–coformer mixture relative to the pure components reflects the tendency of those two compounds to cocrystallize. Thus, predictive calculations may be performed with decent effort on a large set of molecular data in order to identify potentially new cocrystal systems. In addition, it is demonstrated that COSMO-RS theory allows reasonable ranking of coformers for API solubility improvement. As a result, experiments may be focused on those coformers, which have an increased probability of cocrystallization, leading to the largest improvement of the API solubility. In a similar way as potential coformers are identified for cocrystallization, solvents that do not tend to form solvates may be determined based on the highest Hexs with the API. The approach was successfully tested on tyrosine kinase inhibitor axitinib, which has a propensity to form relatively stable solvated structures with the majority of common solvents, as well as on thiophanate-methyl and thiophanate-ethyl benzimidazole fungicides, which form channel solvates. © 2012 Wiley Periodicals, Inc. and the American Pharmacists Association J Pharm Sci


INTRODUCTION
Cocrystals can be defined as homogeneous solid phases containing two or more neutral compounds in a crystal lattice with defined stoichiometry, which are solids in their pure form at ambient conditions. 1 The transformation of active pharmaceutical ingredients (APIs) from their pure crystalline form into cocrystals has experienced increasing interest recently. A cocrystal of the API and an additional compound may show modified properties (such as solubility, dissolution rate, and physical and chemical stability) as com-pared with the pure compounds. 2,3 The possibility to improve the bioavailability 4,5 of the API and to create patentable intellectual property 6,7 constitutes a new and highly attractive route for drug development.
Various experimental methodologies are currently employed for cocrystallization including grinding, 8,9 crystallization from melt, 10 traditional solution crystallization approaches, such as solvent evaporation, 11 cooling, or antisolvent addition, and slurry crystallization. 12 These experimental techniques are typically time consuming and expensive. Therefore, the ability to predict the propensity of different coformers to form a cocrystal with the given API is important.
From general consideration, a likelihood of cocrystal formation is related to the miscibility of API and coformer in the solid state. For a crystalline material, miscibility should be defined by the cocrystal lattice energy. In fact, it was demonstrated previously that rationalization of cocrystal formation in certain cases may be achieved by crystal structure prediction techniques. [13][14][15][16] However, these methods are time consuming and cannot be applied for virtual coformer screening. The effectiveness of the current crystal structure prediction methods quickly decreases with an increase in system complexity (number of molecules in the asymmetric unit and their conformational flexibility). Therefore, the majority of the current computational approaches to virtual coformer screening neglect stabilizing long-order packing contributions to the coformers miscibility.
In practice, rationalization of cocrystal formation is typically based on consideration of only dominant contributions to the miscibility. Intermolecular hydrogen (H)-bonding interaction is the most common focus because of its high strength and directionality. For example, the rational design of cocrystallization is typically performed by a crystal engineering approach, which is based on a hierarchy of H-bonded supramolecular synthons. [17][18][19] A theoretical model was recently proposed for virtual screening based on H-bonding propensities of cocrystal formers derived from molecular electrostatic potential surface calculations. 20 Alternatively, statistical analysis of Cambridge Structural Database (CSD) 21 was performed and a model of molecular complementary in cocrystals was suggested for virtual coformer screening, which was based predominantly on shape and polarity of cocrystal formers. 22 H-bond donor and acceptor counts showed no obvious statistical relationship. A potential drawback of the model is that it was trained on cocrystal observations in the CSD database, ignoring potential failures in realistic cocrystal screenings. Babu et al. 23 have calculated Hbond energies of amide and N-oxide synthons at the HF/6-31G * level in order to compute cocrystal formation. In addition, Hansen solubility parameters were recently applied to describe miscibility of API and coformer to predict cocrystal formation to guide cocrystal screening. 24 In the current study, we demonstrate how conductor-like screening model for real solvents (COSMO-RS) 25 fluid-phase thermodynamics computations describing miscibility of cocrystal formers in a supercooled liquid (melt) phase can be applied to virtual coformer screening. It is assumed that the supercooled liquid phase mimics the cocrystal solid state, neglecting long-order packing contributions (an amorphous solid state). An extensive testing of the approach on multiple experimental screening observations, including pharmaceutical APIs paracetamol, bicalutamide, itraconazole, nicotinamide, meloxicam, carbamazepine, and indomethacine, is reported. Concerning the predictivity of the approach, it should be taken into account that a negative experimental result of cocrystallization of an API with a coformer does not completely exclude the possibility that such cocrystal exists. There are many reasons why it just may not have been observed in the specific experimental setup.
In a similar way as potential coformers are identified, solvents that do not tend to form solid solvates (sometimes called pseudopolymorphs) with the API may be found by COSMO-RS fluid-phase thermodynamics computations. The main difference between solvates and cocrystals is the physical state of the isolated pure components: if one component is a liquid at room temperature, the crystals are designated as solvates. Although cocrystal formulation can bring advantages by increasing dissolution rate and bioavailability of the API, solvate formation is typically undesirable. Solvates might be subsequently desolvated in a final drying step of the formulation process. In such a situation, the final polymorph could be metastable and undergo solid-solid transition during its shelf life. In addition, residual solvent levels in the API must be compatible with International Conference on Harmonisation guidelines (http://www.ich.org) in case of incomplete desolvation. Therefore, selection of the solvent system for crystallization which has the lowest probability of forming solvates with the API is important. Such solvent systems in general may be used directly for slurry crystallization of the stable form or for desolvation of the solvated forms by reslurry experiments to facilitate solvent-mediated transformation and conversion to a stable anhydrous nonsolvated form.
In the current study, we demonstrated how fluidphase thermodynamics calculations allow selecting solvents, which do not form solvates with a tyrosine kinase inhibitor axitinib (trade name Inlyta), and with fungicides thiophanate-methyl (TM) and thiophanate-ethyl (TE).
Previously, such an idea was suggested 26 but never attempted for the solvent selection using the Fábián 22 model.

Approach
COSMO-RS (COnductor-like Screening MOdel for Real Solvents) is a universal theory to predict the thermodynamic equilibrium properties of liquids, which was originally developed by Andreas Klamt. 25,27 COSMO-RS thermodynamics is based on the statistical physics of interacting molecular surface segments. The polar and H-bond interaction energies are quantified based on the surface screening charge densities, which result from a quantum chemical continuum solvation calculation. Because of its ability to treat mixtures at variable temperatures and to compute accurate solvation energies based on first principles, it has become very useful in chemical engineering and in wide areas of physical and medicinal chemistry.
A complete computational modeling and prediction of the cocrystallization process is currently out of reach due to the complexity of the involved steps such as nucleation and crystal growth. However, with COSMO-RS being a fluid-phase thermodynamics model, it is possible to compute a virtually supercooled liquid mixture of the cocrystallization components and obtain the excess enthalpy (H ex ) of stoichiometric m:n mixtures created out of the pure components A and B: H pure and H AB represent the molar enthalpies in the pure reference state and in the m:n mixture, with mole fractions x m = m/(m + n) and x n = n/(m + n). H ex contains all enthalpic contributions and is not limited to H-bonding interactions, although those may be separated from the overall enthalpy by COSMOtherm software. 28,29 We found that H ex is a superior descriptor to the pure H-bonding interaction. Compounds with H ex < 0 are strongly attractively interacting in solution (supercooled liquid) and prefer the mixture enthalpically over their pure liquids. In an extensive set of test calculations presented below, we demonstrated that coformers miscibility as measured by H ex corresponds nicely with an increased probability of forming cocrystals. Since it is reasonable to assume that enthalpic preference of such supercooled liquid phase will also pertain in a mixed crystal, it is plausible to use the liquid phase H ex as a guide for cocrystal screening.

Methods
Excess enthalpies, H ex , were calculated by the COSMOtherm software. 28,29 The screening charge densities for COSMOtherm calculations were generated by the Turbomole package, 30 using the BP86 density functional 31-33 with a TZVP 34 basis set (BP-TZVP-COSMO level of theory). The COSMOfrag 3.3 module 35 was adopted for increased screening performance of large dataset of coformers. Multiple conformations of APIs and coformers were generated by COSMOconf 2.1 36 or OMEGA 37 (YAA) softwares and adopted for the H ex calculations.

Test of the Approach
For an extensive testing, experimental results of cocrystal screening for multiple APIs were taken from literature sources. 19   Overall performance of the COSMO-RS model for coformer selection was estimated by receiver operator characteristic (ROC) curves (Fig. 1). An ROC curve plots the sensitivity (number of true positive predictions/total number of positive observations) versus 1-specificity (number of false-positive predictions/ total number of negative observations) for a binary classifier system (cocrystal screening results) as its discrimination threshold (H ex cutoff) is varied from small to higher values. The area under the curve (AUC) measures the overall performance of the model. Predictions with higher AUCs are generally better and should always be higher than 0.5, indicating the model is better than random selection.
Sometimes the goal of solid crystalline formulation is to find at least one cocrystal former. That is why in addition to the ROC curves, which measure the overall performance of the model, we used enrichment factor (EF) criteria, which measures the probability of getting the first hit (cocrystal) based on H ex ranking relative to the random selection. EF in this case is defined as follows: EF = (1/N fh )/(n h /N). Here, N fh is the number of coformers screened to get the first hit; n h , Figure 1. ROC curves of virtual coformer screenings based on H ex for (a) 3-cyanophenol, 19 (b) 4-cyanophenol, 19 (c) 4-cyanopyridine, 19 (d) bicalutamide, 19 (e) itraconazole, 38 (f) nicotinamide, 40 (g) indomethacin, 24 (h) indomethacin, 39 (i) benzamide, 42 (j) meloxicam, 41 and (k) paracetamol. 43 Enthalpies H ex are calculated at the BP-TZVP-COSMO level of theory. total number of hits [coformers forming cocrystal(s)] in the screening set; N, total number of coformers in the screening set. The maximum enrichment factor, EF max (N fh = 1), is equal to N/n h .
As an example, we present here details of a rational coformer selection for itraconazole cocrystallization. Results of crystal engineering of cocrystals of antifungal drug itraconazole with 1,4-dicarboxylic acids were recently reported. 38 The cocrystals dissolution behavior was shown to be more similar to commercial Sporanox beads (amorphous) than to micronized crystalline itraconazole. The COSMO-RS selection of the coformers based on the H ex s displayed an excellent performance ( Table 1). All experimentally observed cocrystals are ranked at the top of the list with no false-negative outliers. The EF displayed the maximum value of 2. Please note also that even the difference toward cocrystal formation between the transand cis-isomers (see Discussion), fumaric and maleic acid, respectively, are predicted correctly at the BP-TZVP-COSMO parameterization.
The performance of coformer screening based on the H ex as computed with the COSMOtherm program at the BP-TZVP-COSMO level is presented by ROC curves for all the test cases considered in the current study in Figure 1. The corresponding AUC and EF values are summarized in Table 2.
The overall performance of the proposed model is quite good. In seven out of 12 virtual coformer screenings, the AUC was equal or higher than 0.7. EFs for eight coformer selections displayed the maximum values, EF max . A poor performance was observed for reproduction of indomethacin experimental screening observations. That probably can be accounted for by relatively strong contributions of lattice packing effects in the indomethacin cocrystals, which are ignored by the COSMO-RS calculations. Furthermore, most of the cocrystals observed in Ref. 39 were obtained by transformation from the solvate with dioxane (H ex = −1.3 kcal/mol). It cannot be completely ruled out that using a different, less strongly bound solvent would yield additional cocrystals.

Exploring Limitations of the Approach
The proposed method is based on the miscibility of cocrystal formers in a supercooled liquid phase as measured by H ex , ignoring the crystal packing effects. Therefore, optimal screening of coformers should be expected in case differences in H ex contributions exceed variations due to the packing effects.
The above considerations suggest that prediction of differences in the cocrystal forming ability between isomeric compounds may be a challenging task. The surface screening charge densities of isomers are typically almost identical, leading to very close H ex values with a given cocrystal former. If in such a case the packing in the cocrystalline phases is energetically very different, as it may be the case for isomers having differently directed H-bonding groups, resulting in a different H-bonding pattern, a prediction solely based on the mixing behavior will fail.
An illustration of such a case is the cocrystallization screening of isonicotinamide and nicotinamide with 4-hydroxybenzoic acid, clofibric acid, and diclofenac (Table 3) based on a experimental study recently carried out by Báthori et al. 44 It can be seen that the method is getting an ideal enrichment of coformer selection for nicotinamide cocrystallization. However, the obtained H ex results fail to provide any insights into experimentally observed difference in cocrystal forming ability between nicotinamide and isonicotinamide.

Coformer Ranking for Solubility Improvement
Apart from coformer screening for cocrystallization, the COSMOtherm software suite may be applied to coformer ranking for cocrystal solubility improvement. The underlying equations for such calculations were presented previously and are summarized below.
In case of neutral API and acidic coformer, a cocrystal solubility at a given pH may be calculated according to the following equations (see Supplementary material) 45  dissociation constant of an acidic coformer B; : * is the pseudochemical potential of the cocrystal components in solvent (water) or in the pure supercooled liquid state as defined by Ben-Naim. 47 These pseudochemical potentials are evaluated with the COSMOtherm software. Typically, the cocrystal free energy of fusion is unknown and currently cannot be predicted. Therefore, for the solubility ranking, we neglected this contribution, assuming G fus = 0. In Figure 2, the predicted solubility ranking of carbamazepine cocrystals in water is compared with experimental observations. 45 There is a reasonable agreement with experimental data which would allow, for example, selecting nicotinamide coformer for the largest solubility improvement of the drug. The slope between experimental observations and predictions is close to unity within the estimated standard deviation (F).

High-Throughput Computational Cocrystal Screening from Molecular Libraries
The computation of H ex quantities can be done comparatively easy for a large set of compounds in order to identify possibly new API-coformer pairs. As follows, we test such a screening approach for the well-investigated drugs paracetamol and meloxicam by mingling a set of experimentally known coformers with a large list of US Food and Drug Administration-approved compounds from a database called everything added to food in the United States (EAFUS). 48 To generate a molecule library from the EAFUS data we have used the CAS numbers given there and retrieved corresponding SMILES strings from a novel NIH webservice called Chemical Identifier Resolver beta4. 49 From the original EAFUS data all ionic constituents, polymeric constituents with undefined structure, constituents composed of several compounds have been removed resulting in a set of about 2000 compounds. To increase the screening performance, we used the COSMOfrag module, which allows for the rapid generation of an approximate Fprofile for the molecules just starting from a SMILES string based on comparison with a database of already existing F-profiles. 35 In other words, this approach allows generating instantaneously data of nearly quantum chemical accuracy without performing any such explicit and costly calculation. Thus, a screening is done for several hundred compounds within minutes. Figure 2 shows an enrichment plot (rocking plot) for both drugs that have been screened against this list including known coformers of paracetamol 43,50-53 and meloxicam. 41,54,55 The plot has been generated just by ordering the results according to the H ex of each API-coformer pair.
For both drugs, we find nearly all to our best knowledge known coformers within the first third of the screened list. The only exception is the paracetamol coformer naphthalene, which shows up very late, after about 90%; indeed it is quite surprising that such an unpolar compound forms a cocrystal with paracetamol anyway. Its enthalpic interaction with this drug is rather small and therefore this is a case where crystallization is dominated by efficient packing. Concerning meloxicam, the screening obtained sulfuric acid as the compound with the lowest H ex . Although sulfuric acid itself is not known so far to cocrystallize with this API, a literature research revealed that meloxicam hydrogen sulfate is known as a crystalline salt. 56 There is also a somewhat extended plateau in Figure 3, at about 2%-5% of the sampled compounds, which is mainly due to all kind of different nitrogen heterocyclic compounds such as pyridines, pyrazines and oxazoles. A negative H ex of meloxicam with those compounds is caused by the interactions of the meloxicam alcohol group, with the aromatic nitrogen accord-ing to the contact statistics of the COSMOtherm calculations. Those compounds could be interesting targets as alternatives to the acidic compounds used so far in experimental trials.
Of course such a coarse screening will have to be refined by manual selection of the most promising systems by an experienced crystal engineer, but anyway we believe that this is a promising approach to enrich the standard coformer portfolio of drug developers with novel and perhaps unexpected compounds.

Axitinib
Axitinib is a small molecule tyrosine kinase inhibitor developed by Pfizer (Fig. 4a). This API targets the vascular endothelial growth factor to prevent the growth and proliferation of cancer cells via interruption of tumor angiogenesis (formation of vascular supply tissue). Axitinib has shown to be a polymorphically complex API, with five anhydrous forms and 66 solvated forms known. 57,58 Axitinib has a propensity to form relatively stable solvated structures, as a majority of these solvates were characterized as possessing relatively high temperatures of desolvation (desolvation temperatures significantly higher than the normal boiling point of the corresponding solvent), suggesting strong bonding within the crystal. 26,57 These solvates are thermodynamically stable in their corresponding mother liquor and may resist further solvent-mediated transformation to an anhydrous form. From crystallographic consideration, depending on molecular size and H-bonding features, the solvent can either occupy pockets between axitinib dimers forming pocket solvates or link these dimers together by H-bonding with the pyridine or the pyrazole acceptor and the amide donor. 58 Majority of solvates displaying higher desolvation temperatures than the boiling point appeared to be pocket solvates. 26,58 The H ex calculations were performed for 1:1 liquid mixtures of axitinib with 46 pure solvents, described by Campeta et al. 57 (Table 4). Slurry crystallization in all of these solvents except for heptane resulted in solvate formation. 57 Positive H ex values were predicted for 24 solvents, 23 of which form solid solvates with axitinib. That indicates that low miscibility of axitinib with these 23 solvents in the supercooled liquid phase is counterbalanced by lattice packing (threedimensional order) contributions in the solid state, which are ignored in the current calculations. Nevertheless, as in the case of coformer selection for cocrystallization, we assume that miscibility in supercooled liquid as measured by H ex can be sufficient for ranking solvents propensity to form solid solvates.
Heptane was found to have the highest H ex value among all considered solvents (Table 4). That corresponds to the lowest miscibility with axitinib in the supercooled liquid mixture. Heptane indeed was one of the few solvents that did not solvate with axitinib and would be expected to be partially miscible with most solvents at high temperature to facilitate desolvation. In fact, heptane was adopted for solventmediated desolvation experiments through solvates reslurry at 105 • C. 57 The resulting EF provided by H ex ranking is 46.

Thiophanate-Methyl and Thiophanate-Ethyl
Polymorph and solvate formation studies were recently reported for two fungicidal compounds: TM [dimethyl 4,4 -(o-phenylene)bis(3-thioallophanate)] and TE [diethyl 4,4 -(o-phenylene)bis(3-thioallophanate)] (Figs 4b and 4c). 59,60 Although the two molecules are close analogues and differ only by two Both molecules willingly form polymorphs and solvates: TM has two conformational polymorphs and at least 14 solvates, 59,60 whereas four polymorphs and seven solvates were reported for TE. 60 The solvent molecules occupy channels running through the crystal structures, displaying H-bonding and/or van der Waals interactions with the host molecules. Because of the additional CH 2 groups, TE molecules allow formation of larger channels than those observed in the crystal structures of the TM solvates. As a result, majority of the studied solvents do not engage in H-bonding with the TE molecules in contrast to predominance of the H-bonded TM solvates.
Because of a higher mobility of the solvent molecule and the presence of geometric constraints imposed by the cavity size, channel solvate formation seems to be a challenging test to the application based on the H ex calculations only.
Virtual solvent screening results based on the H ex calculations for liquid mixtures of TM and TE with pure solvents are presented in Tables 4 and 5, respectively. Water was found to have the highest H ex value among all considered solvents, supporting lack of the experimental observation of TM and TE pure hydrates. Toluene solvent was correctly ranked as the second top solvent, which does not form solvates with TM (Table 5). Methanol, ethanol, and 1,2-DCE solvents that do not form TE solvates are ranked at the top of the solvent list with toluene as a false-positive prediction ( Table 6). Strong false-negative predictions were given to dimethylacetamide and dimethyl sulfoxide solvents for which TM and TE solvates, respectively, were not observed.
In spite of the complexity of the test, a reasonably good overall performance of the proposed method in application to TM and TE channel solvates is reflected in the maximum EFs of 17 and 12 and AUC values of 0.63 and 0.67, respectively. Experimentally observed solvents which do not form TE solvates are highlighted in gray. 60 Enthalpies H ex are calculated at the BP-TZVP-COSMO level of theory and are presented in kcal/mol. Experimental stoichiometries were used in the calculations whenever available, otherwise 1:1 mixtures were considered.

CONCLUSIONS
It is demonstrated that COSMO-RS theory as implemented in the COSMOtherm software offers a highly efficient way to preselect cocrystal coformers by computational screening. This is achieved by the calculation of the H ex which may be interpreted as the tendency of the two components to associate in the mixture prior to cocrystallization. Most likely due to its detailed and accurate description of all intermolecular interactions, COSMOtherm appears to be more accurate in coformer ranking than some other specially developed procedures, which are focused on intermolecular H-bonding. In addition, it is demonstrated that COSMOtherm allows reasonable ranking of coformers for API solubility improvement. As a result, experiments may be focused on those coformers, which have an increased probability of cocrystallization, leading to the largest improvement of the API solubility.
In a similar way as potential coformers are identified for cocrystallization, solvents which have the lowest probability to form solid solvates may be determined based on the highest H ex values with the API. Such solvent systems may be used directly for slurry crystallization of the stable form as well as for desolvation of the solvated forms by reslurry experiments to facilitate solvent-mediated transformation and conversion to a stable anhydrous nonsolvated form.