Minimizing Polymorphic Risk through Cooperative Computational and Experimental Exploration.

We combine state-of-the-art computational crystal structure prediction (CSP) techniques with a wide range of experimental crystallization methods to understand and explore crystal structure in pharmaceuticals and minimize the risk of unanticipated late-appearing polymorphs. Initially, we demonstrate the power of CSP to rationalize the difficulty in obtaining polymorphs of the well-known pharmaceutical isoniazid and show that CSP provides the structure of the recently obtained, but unsolved, Form III of this drug despite there being only a single resolved form for almost 70 years. More dramatically, our blind CSP study predicts a significant risk of polymorphism for the related iproniazid. Employing a wide variety of experimental techniques, including high-pressure experiments, we experimentally obtained the first three known nonsolvated crystal forms of iproniazid, all of which were successfully predicted in the CSP procedure. We demonstrate the power of CSP methods and free energy calculations to rationalize the observed elusiveness of the third form of iproniazid, the success of high-pressure experiments in obtaining it, and the ability of our synergistic computational-experimental approach to "de-risk" solid form landscapes.

Introduction

Significance of Late-Appearing Polymorphism

Polymorphism is the existence of multiple crystal structures with identical chemical compositions for a particular chemical compound.[1−3] Many pharmaceutical drugs display polymorphism and the different structures have different physicochemical properties such as solubility, hydration stability, etc., which can markedly impact the overall drug efficacy.[4] In addition, factors such as crystal morphology and particle size can influence the drug’s formulation and processing properties.[5] Thus, a thorough understanding of the solid form landscape of a particular compound provides the opportunity to fine-tune material properties.

Polymorph screening, using a wide range of experimental parameters and procedures, has become routine across the solid-state sciences. Famously, McCrone stated that “the number of forms known for a given compound is proportional to the time and energy spent in research on that compound.”[6] Indeed, there are many compounds for which an extraordinary number of forms have been discovered. For instance, the tenth form of galunisertib,[7] the twelfth form of aripiprazole,[8] and the fourteenth polymorph of ROY[9,10] have all been reported. Such a high degree of polymorphism is uncommon, however, and some compounds are monomorphic, i.e., they have only one known pure crystal structure, e.g., Pigment Yellow 74,[11] fenamic acid,[12] and the bromo derivative of ROY.[13] In the ideal case, this is achieved through close packing to give a dense crystal with all potential directional intermolecular interactions being optimal, with no alternatives of comparable stability.[14]

However, it can be premature to assume that a compound is monomorphic based on empirical evidence alone, as there exist a number of examples of late-appearing polymorphs of compounds previously considered to have only a single form, sometimes with significant consequences. Until very recently, isoniazid (ISN) had only one resolved solid form, and the difficulty in obtaining any other crystal structures, despite the molecule being known for almost 70 years, led to it being characterized as monomorphic.[15] Recently, two new forms were discovered viacrystallization from the melt,[16] finally revealing its polymorphism and confirming hints of additional forms suggested by thermomicroscopy experiments decades prior.[17]

The most famous case of late-appearing polymorphism is that of ritonavir, an antiretroviral drug synthesized by Abbott Laboratories with a late-appearing polymorph that resulted in a two-year halt in production and $250 million in lost sales.[1,18,19] Bučar et al. have discussed 10 other cases of elusive or disappearing polymorphs.[1] For example, the dopamine agonist rotigotine, first produced in 1985, was only known to exist in one crystal structure.[20] However, almost 30 years after its discovery, the appearance of a new crystalline form in the Neupro patches stopped its clinical use. The patches were eventually reformulated as a stable amorphous matrix, but during this process, patches were unavailable for four years. Clearly, it is extremely desirable to avoid these late-appearing forms and to be confident that all likely polymorphs of a compound have been discovered and, ideally, fully characterized.

The late appearance of thermodynamic forms may be due to unfavorable crystallization kinetics, but once an energetically preferred form crystallizes, it can inhibit the formation of previously known metastable forms.[21] Interestingly, both ritonavir and rotigotine are examples of conformational polymorphism, in which different molecular conformations are adopted in the polymorphs and the need for conformational change may contribute to the nucleation barrier for a particular form.

The discovery of novel forms is aided by increasingly sophisticated crystallization techniques, e.g., the templating of the catemeric Form V of carbamazepine via sublimation onto the isostructural dihydrocarbamazepine,[22] the crystallization of novel β-coronene under an external magnetic field,[23] and the crystallization of two novel forms of adefovir dipivoxil in ionic liquids.[24] These kinds of experiments are unlikely to feature in a traditional polymorph screen, and hence, combined experimental and computational modeling approaches using a broad range of techniques are needed to minimize the risk of a late-occurring, stable solid form.

It is in this context that computational methods for crystal structure prediction (CSP) can be employed to understand these risks.[25−28] At a basic level, predicting, enumerating, and ranking the stable crystalline forms of a given molecule can provide a picture of what crystal forms are possible and whether the most highly ranked structures (by a chosen scoring function, typically the lattice energy) have all already been experimentally characterized. The results of CSP can, therefore, motivate additional effort in exploring crystallization conditions.[29] At a more sophisticated level, there is the possibility of guiding the experimental polymorph search, both by predicting new forms’ existence and by suggesting conditions under which they might be produced, either to streamline efforts to obtain them or to highlight conditions that should be avoided to minimize the risk of their formation.

Computational approaches can also offer insight into the likelihood of the formation of novel polymorphs under nonambient conditions, particularly high pressure, as was performed post hoc for 2-fluorophenol and 4-chlorophenol,[30] and more recently used a priori to predict high-pressure polymorphs of the pharmaceutical dalcetrapib.[21] However, experimental crystallization of high-pressure polymorphs is often not thermodynamically controlled,[31] frustrating direct comparison between computed high-pressure thermodynamics and high-pressure crystallization experiments. Calculation of crystal structure free energy rather than static lattice energy[32−35] can close the gap between the simulation environment and experimental conditions, but these require expensive and sensitive dynamical calculations.

It is well-known that the “energy landscapes” provided by CSP, the set of predicted stable crystal structures, ranked by their static lattice energy, feature many more unique structures than have been observed even for highly polymorphic molecules like ROY.[36] Clearly, not every structure predicted by a CSP procedure is a feasible polymorph, so to be useful, CSP must suggest which structures are in the “danger zone” on the landscape, i.e., those that are likely to crystallize alongside or instead of the known or desired form(s).

Rationalizing Polymorph Risk through Computed Energetics

To define this “danger zone”, previous work has demonstrated that computed lattice energy differences between experimentally observed polymorphs for a wide range of organic molecules are usually smaller than 2 kJ/mol, less than 7.2 kJ/mol in 95% of cases, and extending beyond 10 kJ/mol in only rare cases.[32] These statistics can be applied to assess the likelihood of observing predicted polymorphs of a target molecule based on their lattice energy difference relative to the lowest-energy structure. From an experimental perspective, this energy window could be measured relative to the lowest energy observed structure. However, this makes the interpretation of CSP results dependent on the extent of experimental screening and thus system-dependent and subject to change if a lower energy structure is observed. With the results of CSP in hand, assuming that the computational exploration for structures is complete and that the energetic ranking is accurate, we know the structure and the energy of the lowest energy possible crystal structure. This global energy minimum should always be assumed to be an observable crystal structure, even if the kinetics of crystallization to this structure are unfavorable. Thus, where CSP has been performed, the most consistent choice of reference for the energetic window of polymorphism is the predicted global minimum, ensuring that the analysis is applied consistently where CSP has been performed “blind”, in advance of any known structures, and to systems with crystal structures already known.

While the main focus of polymorph risk analysis is normally the identification of polymorphs that are thermodynamically more stable than structures that are already known,[7,37] higher energy polymorphs are of interest for a complete understanding of a molecule’s solid form diversity. Desolvation of solvates has been demonstrated as a route to high energy polymorphs, but their instability and difficulties in obtaining high quality crystals of these forms may make them underrepresented in studies of the energy range of polymorphism, which rely on fully determined crystal structures. The relative energies of such desolvated structures can vary widely. Forms II and III of galunisertib, for example, are 8–10 kJ/mol higher in energy than the (unobserved) CSP global minimum and only accessible via desolvation.[7] In other cases, desolvated structures’ relative energies can range from +15 kJ/mol[38] to +25 kJ/mol,[39] and even up to +50 kJ/mol in the case of porous organic frameworks.[40] Clearly, if desolvation or other high energy processes for obtaining alternative forms are under consideration then the energy window considered relevant in CSP must widen significantly.

Regardless, to make any confident prediction of monomorphism (or completely characterized polymorphism), the sampling of possible structures must be sufficiently extensive to have found all the relevant (low-energy) candidates and the method for obtaining their relative energies must be sufficiently accurate.[25] The sampling problem is one of the biggest challenges of CSP due to the high dimensionality of the search space and is made even more difficult when molecular flexibility adds to the number of degrees of freedom. A routine, rudimentary approach for treating flexibility is to sample the crystal packing possibilities of multiple molecular conformers generated from isolated-molecule optimizations. This approach assumes that the in-crystal conformation is nearly equivalent to a stable gas-phase conformer but allows for less stable conformers that might lead to more favorable packing interactions to be considered in the CSP procedure. The advantage of this CSP approach is that the cost of the procedure increases only linearly with the number of conformers considered; the disadvantage is that it does not allow for significant conformational distortions of a gas-phase minimum that might be stabilized (or kinetically trapped) by subsequently favorable packing arrangements.

The severity of this rigid-molecule approximation can be relieved by optimizing the final crystal structures using a method that allows intramolecular degrees of freedom to relax. The most commonly employed such method is periodic density functional theory (DFT).[25] Refinement of predicted structures with DFT also often provides improved energetic rankings and computed properties of the structures compared to the force fields used in the initial structure generation and minimization, albeit at a considerably increased computational cost.

In this work, we combine CSP with both traditional and nontraditional experimental crystallization approaches to explore the polymorphism of two related molecules, ISN and iproniazid (IPN), Figure . As described, until very recently,[16] only one nonsolvated crystal structure of ISN was reported despite screening; our previous, more targeted attempts at gel-assisted crystallizations (with two mimetic gelators) as well as microemulsion crystallization experiments did not obtain any new forms.[15] Even so, there were hints of other possible forms (from the aforementioned thermal microscopy experiments) that had not been further described in almost 50 years.[17] The recent work of Zhang et al.[16] has located two metastable forms (Forms II and III) of ISN via crystallization from the melt. Of these, only Form II’s structure could be fully solved, and a unit cell was proposed for Form III from powder X-ray diffraction data. The elusiveness of these forms of ISN in previous experiments, combined with the tractability of ISN from a CSP perspective as a small, conformationally rigid molecule, make it an ideal candidate for further experimental screening combined with computational study.

Figure 1

Molecular structures of isoniazid (ISN) and iproniazid (IPN).

The structural analogue IPN was the first monoamine oxidase inhibitor and the isopropyl group adds to its conformational flexibility. This may give rise to conformational polymorphism in the solid state in a manner that is unavailable to ISN, as well as possible differences in observed hydrogen bonding motifs due to having one fewer terminal hydrazine hydrogen compared to ISN. Thus, we were interested in how this chemical change impacts the crystal structure landscape. Furthermore, IPN has no reported crystal structures beyond a phosphate salt[41] and is therefore an excellent candidate with which to undertake a fresh cooperative experimental and computational polymorph screen more akin to the screening process in the pharmaceutical industry.

We aim to predict and characterize the polymorph landscapes for these compounds as completely as possible. We begin by predicting “blindly” (i.e., with no prior crystal structure information) the CSP landscape of ISN and IPN. Armed with this knowledge, we attempt to crystallize as many as possible of the forms that appear from the CSP landscapes to be experimentally “at risk” of formation. We determine this risk based in part on the 7.2 kJ/mol energy window for likely polymorphism but with particular emphasis on locating the global energy minimum predicted structure. Our experimental screen incorporates a wide range of techniques, typical solvent screening methods, templated sublimations, gel-phase crystallizations, and high-pressure experiments, to maximize our ability to locate any elusive, metastable, or previously unobserved polymorphs. Finally, we combine the information from both approaches to characterize the observed forms and assess the risk of late-appearing polymorphism for both compounds.

Computational Methods and CSP Procedure

Generation of Hypothetical Structures via CSP

Initial molecular conformers for both molecules were generated via a combined molecular mechanics sampling and DFT optimization procedure (described more fully in the Supporting Information). This procedure yielded a single conformer for ISN, and for IPN a set of five conformers labeled A through E in order of decreasing relative stability (see the Supporting Information for diagrams). All five conformers lay within 5 kJ/mol of the global gas-phase minimum, i.e., with conformational energies well within the expected energy bounds for observed crystal structures.[42]

For each conformer of each molecule, hypothetical crystal packing arrangements were generated with rigid molecular geometries, in our global lattice energy explorer method,[43] details of which are described fully in previous work. The search was restricted to the 25 most common space groups (see the Supporting Information) observed in the Cambridge Structural Database[44] for one molecule in the asymmetric unit (Z′ = 1) and, for isoniazid, also Z′ = 2. Such restrictions on the complexity of the asymmetric unit are a crucial assumption in CSP for maintaining tractable computational cost; however, it precludes the identification of higher Z′ structures and those in unusual space groups. The recently discovered ISN Form II has unusually low symmetry, with four independent molecules in the asymmetric unit (Z′ = 4) and represents a considerable challenge for any CSP effort, due to the number of independent degrees of freedom. Given that there was limited evidence of any polymorphs of ISN at the outset of our work, a search up to Z′ = 2 represents a reasonable compromise between exploration and affordability.

The hypothetical packing arrangements were optimized in a multistage process of successively higher-accuracy energy minimization methods, progressing from rigid-molecule pairwise force field models using the DMACRYS software[45] (see the Supporting Information for a complete description) to a periodic DFT optimization and ranking of the most stable structures. Duplicate structures were removed by automated comparison of computed PXRD patterns obtained via the PLATON[46] program. For the final optimization in periodic DFT, we used the PBE functional[47] and Grimme’s D3 dispersion correction[48] with Becke-Johnson damping (GD3BJ)[49] and a plane-wave basis set, using the VASP[50−53] software package. All atomic degrees of freedom and unit cell parameters were relaxed in the final stage of the procedure, which introduces a description of the molecular flexibility in response to the crystal packing arrangement.

Predicting Free Energies

The energy landscapes obtained via the above methods are computed under the assumption of a static crystal structure; no thermal effects or zero-point energy are included. This is a significant approximation but necessary to manage computational cost. Once a sufficiently small number of plausible structures have been identified, it becomes tractable to predict free energy differences. For ordered crystal structures, the most important contribution to free energies beyond the static energy arises from the dynamics of the crystal structure: the zero-point vibrational energy and the phonon modes populated at finite temperature, which contribute the vibrational term Fvib(T) to the Helmholtz free energy A(T):Fvib(T) is calculated from the phonon frequencies derived from DFT (PBE-GD3BJ). We employed the Phonopy[54] package to obtain the phonon frequencies via finite displacements, calculating energies in VASP for a supercell of each structure; the supercell dimensions are chosen to correspond to sampling reciprocal space q-points of at most 0.12 Å–1 spacing, which has been demonstrated[32,33] to be sufficiently converged for polymorph vibrational energy differences. Phonopy employs a variant of the Parlinski-Li-Kawazoe method[55] for interpolating between explicitly sampled q-points. It is critical to ensure that structures are at true minima to obtain reliable frequencies of vibration about the atomic equilibrium positions; hence, structures were reoptimized with significantly more stringent convergence criteria (1000 eV basis set cutoff, 0.005 eV/Å in forces) before the dynamical matrix was calculated.

Free energy calculations were improved by employing the Debye model to describe the acoustic mode contribution to the phonon density of states from the Brillouin zone center to the nearest sampled q-point.[33] We obtain the elastic tensor from the DFT calculations and calculate an orientationally averaged velocity of sound in the crystal, from which a Debye frequency was determined and used to calculate a correction term for long-wavelength acoustic contributions to Fvib(T). Details are provided in our previous work.[33]

Experimental Methods

Experimental methods were undertaken continuously in unison with the computational studies.

Solution Crystallizations and Gel Phase Crystallizations

Solution crystallizations were performed by the heating of a saturated solution of either ISN or IPN until completely dissolved. The solutions were left to cool slowly in a heating block. These were carried out in parallel with gel-phase crystallizations under the same conditions but in which the heated solution was used to dissolve the gelator (1 w/v%). Then the solutions were also left to cool slowly in the heating blocks.

Templated Sublimations

The powder of the sample being sublimed was placed on a glass slide and then on a Linkam LTS420 heating stage. The templating crystal was affixed to a borosilicate glass coverslip with a small amount of Vaseline and then separated from the glass slide with a small rubber O-ring. The powders were then heated to a temperature to achieve sublimation, 131 °C for IPN and 141 °C for ISN at either 5 or 10 °C/min for 6 h to ensure that all the sample sublimes and then the sample was left at the set temperature overnight to allow for crystal growth.

High-Pressure Experiments

High-pressure experiments were conducted by compressing crystals that were grown at ambient pressure in a Merrill–Bassett diamond anvil cell (DAC)[56] using Fluorinert FC-70 as an inert pressure transmitting fluid. A 250 μm thick stainless steel gasket was preindented to ca. 150 μm and drilled with a 300 μm precision hole to create the sample chamber between the two diamond anvils, culet size of 800 μm. The pressure inside the cell was measured after equilibration using a ruby sphere included in the sample chamber by the R1 ruby fluorescence method.[57] The diamond anvil cell was directly attached to a goniometer head and mounted on the diffractometer. Data were collected using the XIPHOS II diffractometer at Newcastle University, using a four-circle Huber Eulerian goniometer with offset chi cradle fitted with a Bruker APEX II CCD area detector and an Ag Kα IμS generator. Data collections of crystals in DACs are poor at locating H atom positions due to shading by the gasket and DAC, which reduces the completeness of the data set.[58−60]

Novel crystal structures of iproniazid are reported; Forms I and II were produced by slow cooling, and Form III was produced by high-pressure experiments. A more detailed analysis of these novel forms, as well as details on gelator synthesis, characterization methods, and computational processes (molecular conformer generation, space group selection for structure generation and structure minimization) and further results of the free energy calculations are all available in the Supporting Information.

Results and Discussion

Isoniazid Crystal Structure Prediction

The CSP results for isoniazid are shown in Figure . The pronounced global minimum energy structure accurately reproduces the historically observed, thermodynamically most stable structure Form I, known since 1954. All of the other structures generated by CSP lie at least 6.5 kJ/mol higher, at the upper reaches of the expected energy window for polymorphism. The global minimum of the Z′ = 2 search is also isostructural to Form I, but in a lower-symmetry space group; relaxing the symmetry constraints has not yielded any lower-energy structures other than Form I. Duplicate structures occurring within the Z′ sets have been removed, but duplicates between sets (e.g., Form I, the global minima in each) have been retained to emphasize that removing the symmetry constraints still locates many of the same low-energy structures as in the Z′ = 1 search, providing reassurance of sufficiently thorough sampling.

Figure 2

CSP landscape for isoniazid (ISN), in which the lowest-energy predicted structure (circled in solid black) matches the longest-known experimentally observed form. Orange crosses indicate structures generated assuming only one isoniazid molecule in the asymmetric unit (Z′ = 1), while blue triangles indicate Z′ = 2. The global energy minimum structure was located in searches with Z′ = 1 and Z′ = 2 (the two points at −134 kJ/mol, 1.48 g/cm3). The pink circle shows where the recently discovered Z′ = 4 Form II lies when optimized using the same methods, while the broken circle indicates the CSP structure whose lattice parameters and computed PXRD pattern closely match those of the unsolved Form III. Energies and densities are obtained from PBE+GD3BJ periodic DFT.

Such a clear thermodynamic preference for an experimentally observed structure in a CSP landscape is unusual,[61,62] molecules typically exhibit multiple low-lying structures on the CSP landscape well inside the energy window of risk (even if these structures have not been observed), and in cases where a crystal structure is known a priori, it is not necessarily the global lattice energy minimum (which may indicate a risk of the spontaneous formation of this more stable polymorph). The existence of Form I as the global minimum in lattice energy and the size of the energy gap between it and other structures (at this level of theory and within the space groups and Z′ values considered) helps to rationalize the previously long-held empirical conclusion that this molecule had only one accessible crystal structure. If ISN were a novel compound undergoing screening, we would suggest on this basis that the risk of late-appearing, stable polymorphs would be low, particularly if high-energy processes such as desolvation of solvatomorphs are not being considered.

The recent work of Zhang et al.[16] in obtaining two new isoniazid polymorphs (Forms II and III) occurred contemporaneously with the present study and refutes these previous conclusions of monomorphism. Form III is a highly metastable polymorph that converts rapidly to Form I and proved too short-lived for full characterization. The unit cell proposed for Form III (3.931(2) Å, 9.754(5) Å, 8.568(4) Å) from synchrotron powder diffraction is in excellent agreement with those of the second-lowest energy predicted structure (a = 3.751 Å, b = 17.164 Å, c = 9.692 Å, β= 95.67, P21/c, Z′ = 1), i.e., the most plausible polymorphic candidate from CSP (Figure , broken circle, 6.5 kJ/mol above the global minimum), apart from cell doubling in one direction (b from our CSP structure is double c from PXRD indexing).[16] The computed PXRD pattern for our CSP structure is also in good agreement with the experimental pattern (Figure ); reflections present in the experimental pattern but absent from our prediction (noted by red arrows in Figure ) correspond to Form I reflections, due to conversion during sample preparation (see the Supporting Information of ref (16)). Hence, we propose that our predicted structure corresponds to the experimental Form III.

Figure 3

Comparison of experimental PXRD pattern (black, originally presented by Zheng et al.) for Form III of ISN and the computed PXRD pattern (blue) for the second-lowest energy CSP structure. Intensities of each pattern have been scaled to the value of the largest peak in each case. The red arrows highlight peaks in the experimental pattern that we do not predict but which correspond to reflections in the PXRD pattern of Form I of ISN.

Form II, with Z′ = 4, would be unlikely to be found in a typical blind CSP attempt. Lattice energy minimization of the structure reported by Zhang et al. places Form II 5.5 kJ/mol above the global minimum (Figure , pink circle), within the expected energy range for polymorphism. The high calculated energies of Forms II and III, relative to Form I, agree with the experimentally observed metastability of these polymorphs, which transform rapidly to Form I, and the energetic ordering agrees with the room temperature stability ranking (I > II > III) from the experiment.[16]

Iproniazid Crystal Structure Prediction

The CSP results for IPN, with Z′ = 1 structures generated from 5 distinct molecular conformers, are shown in Figure . As is common in blind CSP, we posit that the global minimum energy structure (in space group P21/c), which is also the global density maximum, is the most likely to be observed experimentally.[26] A close-lying predicted structure (space group P21) is approximately 1 kJ/mol higher in energy, and we propose as a likely polymorph. These are followed by a gap, then a more varied, nearly continuous distribution of possible structures starting approximately 5 kJ/mol above the global minimum.

Figure 4

CSP landscape for iproniazid (IPN). Data point color denotes space group number, while the shapes of points indicate the gas-phase conformer used to generate the initial crystal structure, with A being the gas-phase global minimum conformation. In subsequent crystallization experiments, Forms I and II of IPN matched the structures in the solid and dashed black circles, respectively. The global static lattice energy minimum, in the dotted circle, would eventually be located via high-pressure experiments as Form III. Energies and densities are obtained from PBE+GD3BJ periodic DFT.

All but 2 of the 13 structures within the characteristic polymorphic energy window of 7.2 kJ/mol feature either the gas-phase minimum conformer A or the fourth-lowest energy conformer D, indicating that the latter can compensate for its higher conformational energy with improved intermolecular interactions or denser packing. Indeed, the global minimum energy (and maximum density) crystal structure features conformer D. Unlike in the case of ISN, there is an alternative structure (featuring conformer A) within a very small energy gap of the global minimum. While not every static lattice energy minimum on a CSP landscape is guaranteed to be an observable polymorph,[36] we consider the presence of two predicted structures based on different molecular conformers at the bottom of the energy landscape as an indication of multiple likely, stable polymorphs of IPN.

Comparison of Predicted Hydrogen Bond Motifs

Isoniazid

The global minimum structure of ISN from the CSP landscape, which matches the experimental structure Form I, features hydrogen bonding between the pyridyl group of one ISN molecule and the terminal amine group of another, forming a C11(8) chain (Figure a). These chains are in turn linked to each other with N–H···N bonds between the terminal amine and the hydrazide nitrogen atom in an adjacent molecule. However, the carbonyl group does not participate in any H-bonding in Form I.

Figure 5

Hydrogen bonding for (a) ISN-I, (b) ISN-II, (c) ISN-III, and (d) the global minimum energy structure of IPN. Blue lines indicate intermolecular hydrogen bond interactions.

The CSP structure matching Form III, in contrast, features the same terminal amine N–H···N(pyridyl) chains instead linked by (hydrazide)N–H···O(carbonyl) bonds (Figure c). While this and several of the predicted higher-energy structures feature H-bonding motifs that involve the carbonyl group as an acceptor (always in a carbonyl-hydrazide H-bond, sometimes sharing the oxygen atom with a carbonyl-amine H-bond), it appears that involving the carbonyl is unnecessary to form a particularly stable crystal structure of ISN, in opposition to what might be expected based on Etter’s rules.[63]

Form II of ISN displays noticeably different hydrogen bonding to Form I or Form III in a more complex arrangement because of its low symmetry (Z′ = 4). Only 2 of the 4 symmetry-inequivalent molecules exhibit pyridyl-acceptor H-bonding, both via a neighboring hydrazide N–H donor rather than pyridyl-amine chains. However, in all 4 molecules the H-bonding involves all the available protons (terminal amine and hydrazide) and the carbonyl oxygen (Figure b).

Iproniazid

In contrast with ISN, many of the low-energy predicted structures (e.g., 12 of 13 structures within the 7.2 kJ/mol window) display hydrogen bonding involving the carbonyl group (most commonly to the hydrazide N–H), suggesting that the carbonyl acting as an acceptor is energetically optimal. The variety of H-bonding patterns available in predicted structures of IPN is reduced compared to ISN by the lack of the accessible terminal amine group, e.g., the pyridyl group appears less likely to participate in H-bonding due to the isopropyl group precluding the formation of end-to-end molecular chains. Instead the second N–H typically H-bonds to other infinite H-bonding chains of IPN (Figure d).

Experimental Crystallization of ISN and IPN

The crystallization of ISN was carried out in 26 solvents via slow cooling of saturated solutions and through the slow evaporation of ISN solutions. In all cases, the known form of ISN was produced, as first reported in 1954.[64] As a result, sublimation experiments were also undertaken as they have previously been used to crystallize metastable forms.[65 −67] ISN was sublimed by heating on a microscope stage below its melting point (171.4 °C) at relatively high heating rates (5 and 10 °C/min). A borosilicate glass coverslip was placed above the sample and the ISN vapor crystallized on the colder coverslip. However, no new forms were obtained.

For iproniazid, slow cooling and slow evaporation in 22 solvents resulted in 14 samples exhibiting diffraction quality crystals. The crystals are all colorless and plate-like (Figure ). Single crystal X-ray diffraction shows that almost all of the crystals are of the same, new form, designated Form I (IPN-I), in the monoclinic space group P21.

Figure 6

Photos of IPN crystals (a) Form I grown from slow cooling a saturated nitrobenzene solution and (b) Form II produced via the sublimation of IPN powder on a borosilicate coverslip. Unit cell (c) Form I and (d) Form II.

In one instance, crystallization from slow cooling of a saturated toluene solution of IPN produced a second form of IPN, designated Form II (IPN-II), in space group Pbca; as with IPN-I, the structure has Z′ = 1. Despite multiple attempts, it was not possible to reproduce this form through solvent crystallization methods, suggesting that IPN-II is likely to have a high barrier to nucleation in solution and is metastable with respect to or is outgrown by IPN-I.[68] Indeed, slurry experiments, in which both forms were added to a saturated solution, resulted in only IPN-I, as confirmed by PXRD (Figure S1), demonstrating that IPN-I is the more stable form under ambient conditions. IPN-II was, however, reproducibly obtained by the sublimation of IPN powder onto borosilicate coverslips.

Comparison of the structural information and geometric overlays using Mercury[69] (Figure S2) of both IPN Forms I and II with the CSP structures reveals excellent matches with two predicted structures, located 5.8 (Form II) and 1.1 (Form I) kJ/mol above the global minimum. Therefore, both these forms of IPN were correctly predicted by CSP, with energies consistent with the statistics for observed polymorphs. The metastable IPN-II is located higher in the energy landscape, while IPN-I matches the second-lowest energy predicted structure; thus, the predicted energetic ordering of these two forms is consistent with the experimental observation that IPN-I is more readily obtained and more stable than IPN-II.

However, these solution and sublimation crystallization experiments did not yield crystals corresponding to the global energy minimum on the CSP landscape. Furthermore, the landscape (Figure ) contains multiple structures that are of very similar energy to IPN-II. It is not apparent why IPN-II is experimentally observed while similarly metastable structures and the global minimum are not obtained. Both observations indicate that the risk of late-appearing polymorphism of IPN remains after conventional solvent screening and sublimation.

Gel-Phase Crystallizations

Gel phase crystallization is an emerging technique for expanding the polymorphism search space[70,71] and was undertaken for ISN and IPN with a series of gelators (Figure ). Gelator 1 mimics the structure of ISN and IPN and parallels the use of other drug-mimetic gelators that have previously been shown to stabilize metastable forms over the thermodynamically favored polymorph,[70] as well as resulting in the discovery of new solvates.[72] Gelators 2 and 3 were chosen to mimic coformers known to form cocrystals with ISN. The cocrystals of ISN are polymorphic and exhibit different H-bonding motifs with the carbonyl and pyridyl groups both being H-bond acceptors.[73−75] Thus, these gelators may be able to template forms with different H-bonding motifs. Gelator 4 forms gels in aromatic solvents and therefore was chosen in an attempt to recover Form II in toluene.

Figure 7

Chemical structures of gelators used in this study.

These experiments were carried out at the same solute concentration as the slow cooling crystallizations. For each experiment, either ISN or IPN and then the gelator were dissolved with heating and sonication and left to cool under ambient conditions. In all cases, Form I of both ISN and IPN were obtained, demonstrating the strong preference for these forms to crystallize over any other predicted structures.

Previous work has demonstrated that it is possible to template the growth of a particular form and discover new polymorphs by subliming onto different surfaces, e.g., polycrystalline powders,[76] siloxane-coated glass,[65] and crystals with related structure.[32,33] Both Forms I and II of IPN display different hydrogen bonding patterns to the known form of ISN and hence crystals of one compound could be used to template the growth of new forms of the other analogue.

In the CSP landscapes, the lowest energy IPN structure with ISN-like H-bonding is 7.0 kJ/mol above the global minimum. Conversely, on the ISN landscape, the lowest energy structure displaying IPN-like H-bonding (i.e., involving the carbonyl group) is the CSP structure that we propose matches Form III of ISN, 6.5 kJ/mol higher in energy than Form I. While these energy differences are toward the higher end of the energy range of expected polymorphism based on experimental statistics, they remain plausible risks. This assessment is borne out in the case of ISN with the recent experimental production of Form III, but it is also significant for IPN, as the observable Form II is only 1.2 kJ/mol more stable than the aforementioned lowest-energy ISN-like structure.

Sublimation of each compound was attempted using crystals of either ISN or IPN as a template for the other by crystallization directly from the vapor phase. Crystals of both compounds did grow on top of the surface of the parent template (Figure S3), but in both cases the same polymorph was produced (Form I ISN or Form II IPN) as sublimation crystallization in the absence of the template. These experiments further qualify the risk of further metastable forms of ISN/IPN with unusual H-bonding patterns; these calculated forms, along with the global minimum form of IPN, appear to be inaccessible through these methods.

We also attempted to seed melt crystallizations with using either ISN or IPN as a template for the other. However, we did not observe changes to crystal morphology (using polarized optical microscopy) or changes of melting temperatures during the reheating measurements and concluded that seeding was unsuccessful. Hence, we did not pursue such experiments further.

Elusive High-Density Form of Iproniazid

The calculated global energy minimum form of IPN is notably denser than either of the forms experimentally obtained. While it is not unprecedented that the first-discovered or most experimentally accessible structure is not the global energy minimum on a CSP landscape, the lower energy and significantly greater density together single this prediction out for further efforts. From the perspective of confidence in solid form screening, an unobserved CSP global minimum represents a major risk and a priority for experimental work. A higher-density polymorph should become comparatively even more stable under higher pressure, suggesting high pressure crystallization as a means to obtain this form experimentally.

Geometry optimizations under pressure show that the energy difference between IPN-I and the global minimum CSP structure widens from −0.7 at zero applied pressure to −7.6 kJ/mol at P = 2.4 GPa (Figure ). This is expected, as the global minimum is already the densest predicted structure and therefore no “crossover” in ranking is expected with increased pressure. While this trend agrees with physical intuition, it provides no clarity as to why this global minimum structure is not found in conventional crystallization experiments. Evidently, static lattice energy calculations alone are not enough to rationalize the elusiveness of the high-density form.

Figure 8

Computed energy difference, ΔE, between IPN-I and the global minimum in lattice energy (reference value) as a function of pressure. Gray ■ are static lattice energies including the PV contribution; the remaining data are Helmholtz free energies presented at three different temperatures: 0 K (black ○), i.e., only ZPE contributions, 100 K (blue +), and 300 K (red ×). Values were calculated via periodic PBE+GD3BJ; vibrational contributions were corrected with a Debye model for low-frequency modes.

Computed Free Energies of Iproniazid

Given the increased possibility of obtaining the elusive high-density predicted structure of IPN under elevated pressures, a more useful quantity to compute is the free energy change as a function of both temperature and pressure. However, the latter is a computationally expensive proposition, as each step in pressure requires new phonons to be calculated as the crystal structure is compressed. To provide some estimate of the free energy difference between the forms under higher pressure, we make the approximation that the pressure and vibrational effects are independent and purely additive, i.e., we compute vibrational contributions to the free energy only for the ambient pressure structures and add these to the PV contribution. These approximate free energies are shown in Figure .

With dynamical effects incorporated into the calculations, the static CSP global minimum becomes metastable with respect to IPN-I at ambient pressure, being 2.8 kJ/mol higher in free energy at 300 K. In fact, the missing high-density form of IPN is computed to be higher in free energy at ambient pressure than IPN-I across all temperatures; the vibrational zero-point energy (ZPE) difference of 1.8 kJ/mol between the two forms reranks them even before any thermal contributions. It is known that vibrational ZPE can be important in relative polymorph ranking, particularly for some hydrogen-bonded species.[77,78] In these cases, an accurate treatment of lattice dynamics is necessary to obtain an estimate of the ZPE and hence polymorph rankings that are consistent with experimental observations, as appears to be the case for IPN.

Having considered both temperature and pressure, the elusiveness of the high-density form is more easily explained: this structure is metastable at standard temperature and pressure with respect to IPN-I, with a considerably larger energy difference between structures than the static calculations alone indicated. As expected, increased pressure stabilizes this high-density form, with higher temperatures requiring higher applied pressures to make it more favorable than IPN-I. While our approximation of additivity of the vibrational and PV terms compounds the uncertainty in the free energy, the results indicate that at 300 K an applied pressure in excess of 1.0 GPa would be required to make the high-density structure most favorable.

Given this reordering, we should consider whether other predicted structures could similarly become competitive with IPN-I in free energy terms with increasing temperature or pressure. The expense and sensitivity of calculating the periodic DFT phonon frequencies makes it desirable to estimate in which cases reordering is likely without carrying out the full calculation on many structures. Our previous work has demonstrated that the vibrational contribution to polymorphic free energy differences rarely exceeds 2 kJ/mol at room temperature.[32] Given that the third lowest energy predicted structure of IPN is 3.8 kJ/mol higher in static lattice energy than IPN-I, any of the higher energy structures are considered unlikely to become more stable in free energy at room temperature, and so we restrict our full free energy treatment to IPN-I and the high-density global minimum structure. (The expected range in magnitude of this entropic contribution also justifies our decision to consider only static lattice energies for ISN, as no predicted structures are likely to become competitive with Form I at room temperature.)

Assessing the polymorphic risk of the IPN landscape based on calculated free energies, the global lattice energy minimum is now less of a risk than indicated by the static calculations alone. This demonstrates the power of computational work to identify risks in crystallization processes but also highlights the importance of thorough, rigorous application of advanced methods to provide accurate insight. If assessments of risk are made on the basis of “black-box” CSP approaches alone, using static lattice energies (as in Figure ), then this maximal density structure would be perceived as a serious risk as a late appearing, stable polymorph because it is the global minimum on such a landscape. This risk is diminished somewhat when free energy (including quantum vibrational effects) is considered, as dynamical contributions rerank the two structures and suggest IPN-I is in fact most stable. However, the free energy difference between the two remains well within the “danger zone” of polymorphism at ambient pressure, and hence, it is prudent to explore whether high-pressure experiments can indeed obtain this form as the free energy trends in Figure indicate.

High-Pressure Experiments

The global lattice energy minimum on the IPN landscape contains the molecule in a less-stable conformation than that of IPN-I or IPN-II. This conformation may have a sufficiently high nucleation barrier that it cannot be crystallized by solution phase, gel phase, or sublimation screening as its crystallization is kinetically hindered. Similar observations were made for ritonavir and rotigotine whereby the thermodynamically favored forms were conformationally different from the metastable forms and were not discovered for many years.

As the CSP static global minimum structure is significantly denser than IPN-I and IPN-II, it is stabilized by higher pressures, due to the smaller pressure–volume contribution to the free energy; this is shown by the free energy calculations (Figure ). High-pressure experiments are known to be capable of effecting conformational change in crystal structures.[79,80] Thus, both high-pressure recrystallization and compression of crystals grown at ambient pressures were undertaken in Merrill–Bassett diamond anvil cells (DACs).

All attempts to recrystallize IPN from solution at various pressures produced polycrystalline samples. This may be due to the many surfaces on which nucleation can occur inside the DAC, e.g., the ruby spheres, the edge of the tungsten gasket, or the faces of the diamond.[81] As a result, we turned to compression of crystals grown at ambient pressure.

At lower pressures (≤0.3 GPa), only slight reductions in the unit cell dimensions were observed for IPN-II (Table S1). However, the compression is anisotropic, affecting mostly the b axis. At higher pressures (ca. 0.5–0.8 GPa), the crystal breaks perpendicular to the longest axis, indicating that this form is unable to compress further or transform to relieve the stress caused by the elevated pressure.

When IPN-I was compressed up to ca. 2.1 GPa, no changes to the crystal habit could be observed visually. However, single crystal X-ray diffraction confirmed that under hydrostatic pressure, IPN-I (P21) undergoes a single crystal-to-single crystal transformation to produce a new structure, IPN-III (P21/c). Pressure-mediated transformations without the destruction or dissolution of the crystal are rare but can be achieved for molecules with conformational flexibility.[81] Single crystal-to-single crystal transformations have been achieved for other organic molecules through conformational changes, e.g., β-glycine to δ-glycine,[82] glutathione-I to glutathione-II,[58] and di-p-tolyl disulfide α form to β form.[83] In the present case, this transformation results in a conformational change such that the structure closely matches the predicted conformer D. Structural overlay (Figure S4) confirms that this new form corresponds to the predicted global lattice energy minimum structure of IPN.

Upon decompression to ambient pressure, the crystallinity of the sample was significantly reduced because of damage to the crystal from compression, decompression, and X-ray irradiation. However, it proved possible to obtain an ambient pressure structure determination, which while of low precision and data completeness unambiguously identifies that Form III is retained after decompression and removal from the DAC. Comparison of the crystallographic data for the Form III high-pressure data and the data for the crystal recovered from the DAC are available in the ESI (Table S2). This key result indicates that Form III represents a considerable risk as a late-appearing polymorph, if created in an industrial setting under milling conditions, for example.

Crystal Packing Analysis

Differences in the packing of the polymorphs of IPN explain why Form I rather than Form II transforms to Form III and why Form III is the densest polymorph. The three IPN polymorphs exhibit a similar H-bonding motif, a C11(4) chain with the carbonyl group as an acceptor and the donor being the N–H adjacent to the carbonyl of the next molecule (Figure and Figure ). These are the shortest contacts for all forms, and they appear as two spikes on the 2D fingerprint plots of the IPN Hirshfeld surfaces (Figures S5 and S6), which were used to visualize the differences in intermolecular interactions. For IPN-I, the H-bonding chains run parallel to the b-axis and only a small compression of the b-axis is observed upon transition to IPN-III (Table ). This is commensurate with only a small reduction in the H-bond length in IPN-III (Table S2). Similarly, these chains run parallel to the a-axis for IPN-II and only a modest reduction of this axis is recorded after compression in the DAC (Table S1). This is in line with previous studies that the shortest contacts remain unchanged and transformations instead rely on a rearrangement of longer contacts.[84,85]

Figure 9

Hydrogen bonding for (a) IPN-I, (b) IPN-II and (c) IPN-III. Blue lines indicate interactions between molecules.

Figure 10

H-bonding sheets of IPN visualized perpendicular to the H-bonding chains with alternative sheets colored red and blue. Blue lines represent interactions between IPN molecules. (a) IPN-I, (b) IPN-III, and (c) IPN-II.

Table 1

Crystallographic Data for the Three Polymorphs of IPN

crystal form	Form I	Form II	Form III
formula	C9H13N3O	C9H13N3O	C9H13N3O
molecular weight/g mol–1	179.219	179.219	179.219
crystal system	monoclinic	orthorhombic	monoclinic
space group	P21	Pbca	P21/c
T/K	120	120	291
pressure/GPa	ambient	ambient	2.21
a/Å	8.1440(8)	4.9971(7)	11.240(5)
b/Å	5.0966(5)	16.831(3)	5.043(3)
c/Å	11.7207(13)	22.850(3)	15.171(14)
α/°	90	90	90
β/°	107.122(4)	90	109.14(7)
γ/°	90	90	90
V/Å3	464.93(8)	1921.7(5)	812.4(10)
Z	2	8	4
Z′	1	1	1
ρcalc/g cm–3	1.280	1.239	1.465
independent reflections	3067 [Rint = 0.0731]	1577 [Rint = 0.0918]	287 [Rint = 0.0851]
goodness-of-fit	1.024	1.060	1.082
final R indexes [I ≥ 2σ(I)]	R1 = 0.0501	R1 = 0.0398	R1 = 0.0987
	wR2 = 0.1033	wR2 = 0.0857	wR2 = 0.2524
final R indexes [all data]	R1 = 0.0719	R1 = 0.0612	R1 = 0.1434
	wR2 = 0.1111	wR2 = 0.0945	wR2 = 0.2863
CCDC code	2011025	2011026	2011027

Accordingly, there are some subtle differences between forms in the interactions of the pyridyl groups in adjacent H-bonding chains. All IPN forms exhibit a short contact between the pyridyl groups via a (pyridyl)C–H···N(pyridyl) interaction. For IPN-I and IPN-II, these are C11(3) chains, which stack in a herringbone-like manner. These contacts in IPN-II are longer than in IPN-I (Table S3), indicating poorer packing. The compression of IPN-I significantly reduces the a-axis (the c-axis in IPN-III), such that the pyridyl groups are less offset from one another. The conformational change also results in the two pyridyl rings becoming nearly coplanar. These changes transform these contacts to R22(6) rings in IPN-III which are noticeably shorter than in the other forms (Figure S8g, Table S3); further, the H-bonding chains now run antiparallel to one another (Figure ).

Together, the H-bonding chains and the pyridyl–pyridyl contacts produce sheets that connect the IPN molecules (Figure ) in all three forms. These sheets are additionally held together by (hydrazine)N–H···N(hydrazine) H-bond interactions (Figure S9a and S10). The shorter pyridyl–pyridyl contacts in IPN-III allow the sheets to stack closer together, while the antiparallel chains and the conformational change allows rotation of the isopropyl group for denser packing within these sheets (Figure S9d); both effects contribute to IPN-III being the densest. Further analysis of the packing and presentation of the Hirshfeld fingerprint plots can be found in the Supporting Information.

In IPN-II, the b-axis is normal to the plane of these sheets (Figure S10), and it is this axis along which IPN-II compresses before breaking under pressure. It is speculated that these sheets compress together rather than the H-bonding chains to accommodate the pressure increases.

Though all three forms contain these H-bonded sheets, the shape of the sheets in IPN-I and IPN-III is notably different from those in IPN-II. The sheets of both Forms I and III are relatively planar, with gentle undulations. In contrast, those of Form II are more corrugated, having an increased roughness (rugosity) with neighboring sheets interpenetrating more (Figure ). As these sheets in IPN-II are very different to those in IPN-III, there is no obvious route for transformation to IPN-III, in contrast to IPN-I which only needs to undergo a comparatively small change in the sheet structure. Topological rugosity has been linked to slip planes for crystals which affect the mechanical properties and thus the tablet ability of different forms.[86,87] Therefore, it is expected that the different forms of IPN may behave differently during formulation, particularly if tabletting pressures cause the transformation from IPN-I to IPN-III that can persist at ambient pressure. Given the structural differences, however, we do not expect IPN-II to transform due to pressure during the tableting process.

Conclusions

We have demonstrated the power of computational CSP methods to rationalize the difficulty in obtaining polymorphs of isoniazid, due to the absence of any competitive predicted thermodynamic minima (assuming Z′ ≤ 2). Solvent and sublimation-based crystallization methods did not reveal any further forms beyond the long-known Form I, consistent with previous screening.[15]

The contemporaneous discovery by Zhang et al. of two new forms from melt crystallization[16] is remarkable but consistent with our computational findings in both cases. Form II would not be predicted in a typical CSP due to its high Z′, making a search prohibitively expensive in the general case. However, the structure of Form II is ranked favorably in energy when treated with our computational methods, and within the range of likely energy differences between observed polymorphs. Similarly, the second-lowest energy CSP structure of ISN lies somewhat higher in energy but still within the usual range for polymorphism and matches the PXRD pattern of Form III, which could not be fully characterized experimentally. Therefore, although our initial assumptions precluded our CSP procedure from predicting Form II, Form III was predicted as the most likely polymorphic structure, albeit with an energy that suggests it would be challenging to isolate (and arguably agreeing with its experimental instability). Thus, our approach successfully predicted the only other Z′ = 1 structure of ISN yet discovered and has revealed its crystal structure.

More significantly, a blind CSP of the analogue iproniazid predicted a polymorphic system with at least two notably low-energy structures and several metastable ones. An exhaustive experimental screening process, including solvent-based, gel-phase, and sublimation crystallization, successfully obtained two of the predicted structures: the stable Form I and metastable Form II. However, the global minimum on the static energy landscape, the densest predicted structure, remained elusive despite its predicted thermodynamic stability. It was only through diamond anvil compression experiments that this structure, Form III, was obtained experimentally. Detailed free-energy calculations representing the state-of-the-art in CSP techniques rationalize the stability relationship between Forms I and III. Form III is not the global energetic minimum when thermal effects are considered but is obtainable when experimental searching is guided by CSP. Once formed by compression, Form III persists at ambient pressure and hence CSP has revealed a high-risk late appearing polymorph and indicated the conditions under which it is formed.

This work demonstrates the power of combining exhaustive experimental screening with modern CSP methods to elucidate the risk of late-appearing polymorphism. We emphasize the synergy between the two fields; CSP of iproniazid suggested a significant risk of polymorphism, which justified a thorough experimental screening that obtained Forms I and II. When the most likely structure according to CSP eluded this screening, further computational analysis via free energy calculations rationalized its elusiveness, while high-pressure experiments motivated by the predicted maximal density of this structure successfully yielded Form III. The crystallization and characterization of three polymorphs of iproniazid, for which no pure crystal structure was previously available, is tangible evidence of the value of this combined approach for exploring and “de-risking” solid form landscapes.