Review pubs.acs.org/OPRD
Current Computational Approaches to Support Pharmaceutical Solid Form Selection Yuriy A. Abramov Pfizer Global Research and Development, Groton, Connecticut 06340, USA ABSTRACT: This contribution reviews current computational approaches to complement polymorph screening in the pharmaceutical industry by facilitating an assessment of the likelihood of a missed stable form. The focus is given to such methods as crystal structure prediction, hydrogen bonding propensity analyses, and rational solvent selection for stable form screening. Special consideration is given to the advantages and limitations of each method. It is suggested that, due to the outlined limitations, an in silico risk assessment of solid form selection in the pharmaceutical industry should be based on a combination of all of these approaches. A detailed application of these computational tools in the pharmaceutical industry is illustrated on a selective dual ALK and c-MET inhibitor crizotinib (trade name Xalkori) and on a VEGF inhibitor axitinib (trade name Inlyta).
1. INTRODUCTION Polymorphism of the crystalline state of pharmaceutical compounds is quite a common phenomenon, which has been the subject of intense investigation for over 40 years.1 The scope of this phenomenon in the pharmaceutical industry was recently demonstrated by an analysis of a series of 245 polymorph screens performed at SSCI (http://www. ssci-inc.com) indicating that about half of the compounds exhibited polymorphism.2 It was also suggested that approximately one-third of organic compounds and about 80% of marketed pharmaceuticals exhibit polymorphism under experimentally accessible conditions.3 In the pharmaceutical industry, drug polymorphism can be a critical problem and is the subject of various regulatory considerations.4 Polymorphism can impact the performance of a drug in terms of solubility and bioavailability,5 chemical and physical stability,6 and mechanical properties.7 That is why polymorph screening is typically performed to identify thermodynamically stable crystal forms under normal manufacturing and storage conditions to ensure that the manufactured form does not undergo changes during the life cycle of the drug product. The importance of thermodynamically stable form selection in the pharmaceutical industry can be illustrated by well-known examples of polymorph-induced impacts on marketed drugs Norvir (ritonavir) and Neupro (rotigotine patches). In the first case, Abbott Laboratories had to stop sales of Norvir in 1998 due to a failure in a dissolution test, which was caused by the precipitation of a more stable and less soluble form II of the compound.8 In the second example, undesirable crystallization of rotigotine was found in the patches that were used to administer the drug. These crystals formed as a result of a previously unknown stable polymorph of rotigotine, causing UCB to suspend the marketing of this drug in the US.9 An example of process development challenges to accommodate a late-appearing polymorph for a fast-track development compound was also reported by BMS scientists.10 This contribution reviews current computational approaches to complement and direct stable form screening in order to © 2012 American Chemical Society
assess the likelihood of missing a stable form at early stages of drug development. In particular, a focus will be given to such methods as crystal structure prediction (CSP), hydrogen bonding (H-bonding) propensity analysis, risk of conformational polymorphism, and rational solvent selection for stable form screening. Special attention will be paid to limitations of each method. The major question which these computational approaches try to address is as follows: “How likely is an alternative stable form?” in order to understand whether any additional study should be performed or a polymorph screening, which had not discovered a new stable form, was likely to be complete with no additional study being required. Consideration of such important topics as computational support of crystallization and polymorphism of salts and cocrystals,11 solvent, additives, and impurity effects on polymorph morphology12 and prediction of the potential impact on solubility by an unknown stable form5d is out of the scope of this review.
2. CRYSTAL STRUCTURE PREDICTION Most logical and at the same time the most complicated theoretical approach to study drug polymorphism is based on the crystal structure prediction (CSP) approach. CSP is a rapidly growing area of theoretical research and development with a history of more than 20 years. In spite of obvious potential benefits for different industries, including pharmaceutical, the initial goal of CSP was limited to the development of computational approaches to correctly predict an observed crystal structure from the chemical diagram of the compound. A typical assumption was made (usually with no reasonable polymorph screening performed) that this crystal structure corresponded to the most stable form. In addition, the greater focus of CSP until recently was given to relatively simple model systems to allow verification of various assumptions and Special Issue: Polymorphism and Crystallization 2013 Received: October 2, 2012 Published: December 11, 2012 472
dx.doi.org/10.1021/op300274s | Org. Process Res. Dev. 2013, 17, 472−485
Organic Process Research & Development
Review
likely to be complete and no additional study is needed. 3Oxauracil provides an example of this type of lattice energy landscape with the observed crystal structure being more stable than any other computationally predicted structures by 4 kJ/ mol.18 In that case a limited crystallization screen did not reveal any additional form except for a monohydrate. The crystal energy landscape that suggests the possibility of an alternative crystal form, which is more thermodynamically stable than the known form and has yet to be found, is shown in Figure 1c. This pattern is indicative of the potential risk of solid form selection and would suggest additional polymorph screening to find the most stable form. This type of information is very useful in helping to avoid late stage thermodynamically stable form discovery, such as form II of ritonavir or rotigotine. 5-Fluorouracil provides an example of such a type of predicted lattice energy pattern leading to the discovery of a new form that corresponded to the global energy minimum structure from the computational search.19 Finally, the polymorphic crystal energy landscape, displaying many minima close in energy, is presented in Figure 1d. Such a landscape may provide an indication of a strong competition between different intermolecular interactions leading to a difficulty of packing in the way which optimizes all of the interactions. This pattern also indicates a potential risk of solid form selection and would suggest additional polymorph screening to find the most stable form. An example of such a lattice energy landscape is provided by aspirin, which had two stackings of the same sheet almost equi-energetic at the global minimum.20 Despite the visible simplicity of the scenario of CSP outcomes, as presented in Figure 1, an interpretation of crystal energy landscapes should be performed bearing in mind the challenges and limitations discussed below. Potential Energy. The values of lattice energy difference between the computationally predicted structures depend on the level of theory applied. The energy difference between possible structures are typically very small so that minor changes in the way the energy is calculated could lead to a significant reranking of the predicted structures. Therefore, special attention should be given to the selection of the model (level of theory) for lattice energy calculations, which should reflect a balance between accuracy and computational expense. It has been shown in multiple studies of small rigid molecules that an important improvement in the reliability of crystal structure prediction for polar H-bonding molecules is achieved through the use of a more elaborate representation of electrostatic interactions. Among the most successful methods there appears to be the density functional theory with dispersion energy corrections (DFT-D)21 and a model with distributed multipole expansion (DMA) of electrostatic interaction.22 Temperature Effect. Perhaps one of the most important intrinsic limitations of polymorphism prediction by typical CSP techniques is related to the fact that the free energy landscape of possible crystal structures at ambient temperature is approximated by a static lattice energy landscape at 0 K, ignoring zero-point vibrations. Even in the case of the application of force fields parametrized to predict structures at ambient temperature (e.g., Williams 99 force field)23 the entropic contribution to the polymorphic free energy would not be described. However, only entropic contributions allow accounting for an enantiotropic relationship between polymorphs. Therefore, the CSP approach is currently limited to
methodologies. In fact, according to surveys of the CSP literature up to 2004, only 4% of all compounds studied are pharmaceuticals.13 Milestones of method developments and compound complexity progress are reflected in the results of five blind tests14 and are covered in the recent comprehensive reviews.15 This review will focus specifically on consideration of CSP as an important computational tool that is valuable for performing a polymorphic risk assessment for solid form nomination during pharmaceutical development. The currently accepted approaches to interpretation of CSP results for polymorphs are considered. In addition, different limitations of the CSP methodologies, which may affect the interpretation, are discussed. Polymorphism Consideration via Crystal Energy Landscape Consideration. Current CSP techniques describe the possible polymorphism of an organic molecule via lattice energy minimization. Minima of the lattice energy landscape within a certain energy range describe thermodynamically feasible polymorphs of the molecule at 0 K. Barriers of transformation between these minima, which are important for polymorph kinetics, are not considered by the current CSP methods. Price suggested interpretation of different crystal energy landscape patterns in order to study organic molecule polymorphism (Figure 1).13b,16 For example, the crystal energy
Figure 1. Some idealized crystal energy landscapes.13b * is the lattice energy minimum corresponding to the experimental crystal structure. ■ and ● are hypothetical crystal structures, where each ● denotes a distinct crystal structure, with some similarity in the stronger intermolecular interactions to the experimental structure, and each ■ denotes a distinct crystal structure which is so different from the experimental structure that interconversion is likely to be kinetically prohibited. Reproduced from ref 13b. with permission.
landscape presented in Figure 1a demonstrates a crystal form as a global minimum in energy separated from all other structures by a considerable energy gap. The proposed energy gap is about 10 kJ/mol and should depend on the level of theory. Such a crystal energy landscape is an ideal indication that this organic molecule displays only one crystal form and supports a low likelihood of existing as another stable polymorph. Pigment Yellow 74 is one reported example of a true monomorphic case in which the computationally predicted most stable form is separated by an energy gap of 12 kJ/mol from all other simulated structures.17 A more common crystal energy landscape, which corresponds to the most stable form as a global minimum, is shown in Figure 1b. In this case, an energy gap (about or less than 4 kJ/mol) is not considered to be big enough to rule out the possibility of metastable polymorph existence. Such a result from a CSP study would also support the low likelihood of existence of a new stable form and suggests that the polymorph screening, which had not discovered any other stable form, was 473
dx.doi.org/10.1021/op300274s | Org. Process Res. Dev. 2013, 17, 472−485
Organic Process Research & Development
Review
Figure 2. Examples of successful application of the current CSP methods to predict structures of flexible molecules. That includes glycol (a),37 glycerol (b),37 carbamazepine (c),36 phenobarbital (d),38 aspirin (e),20 piracetam (f),39 6-amino-2-phenylsulfonylimino-1,2-dihydropyridine (g),21c,d alanine (h),40 valine (i),40 leucine (j),40 isoleucine (k),40 1,2-dichloro-4,5-dinitrobenzene (l),14e,41 benzyl(4-(4-methyl-5-(p-tolylsulfonyl)-1,3-thiazol2-yl)phenyl)carbamate (m),14e,41 axitinib (n),42 and crizotinib (o).
CSP and metadynamics to study PR179 pigment polymorphism.25 Free energy minimization and metadynamics simulations at ambient conditions effectively identified 7 unstable minima out of 18 low energy structures created by CSP at 0 K. The temperature effect on polymorph relative stabilities moving from 0 K to ambient temperature can be described by two contributions: a difference in heat capacity (which determines the temperature behavior of enthalpy and entropy of each form) and a difference in entropy (−TΔS contribution to the free energy difference of the two forms). Since the difference between heat capacities of polymorphs is normally quite small,26 the difference in polymorph entropy is believed to be a more significant contribution. Entropy is dominated by contributions from low-energy lattice vibrations, which can be modeled by lattice dynamics calculations.26a,27 Vibrational
the study of pharmaceutical polymorphs (monotropically or enantiotropically related), which do not display a solid state transition below or close to ambient temperature. Moreover, even for such systems, thermal motion effects may describe an averaging over multiple lattice energy minima that are separated by low barriers. That would result in fewer minima (i.e., predicted polymorphic forms) in a free energy surface. Therefore, thermal motion effects may, at least partially, account for typically a much lower number of experimentally observed polymorphs relative to the computationally predicted ones. This effect was demonstrated by molecular dynamics simulations of acetic acid at 200 K.24 In that study, many higher energy structures generated by CSP converted into lower energy structures, significantly reducing the number of possible polymorphs. Similar behavior was demonstrated by ZykovaTiman et al. by an application of a combination of standard 474
dx.doi.org/10.1021/op300274s | Org. Process Res. Dev. 2013, 17, 472−485
Organic Process Research & Development
Review
systematic conformational scanning for such a molecule would require a massive and expensive set of calculations. Therefore, a rational selection of the preferred starting molecular conformations for CSP of complex flexible molecules is of great importance to speed up the calculations. Finally, the presence of molecular flexibility introduces an additional contribution beyond rigid body vibrations to the entropic difference between calculated polymorphic forms, which is nontrivial to account for. It was demonstrated that the electronic level calculations are most appropriate for balanced treatment of inter- and intramolecular energy contributions to lattice energies. For this, ideally the lattice energy of a crystal structure would be evaluated entirely quantum mechanically, for example using the DFT-D approach.21a Alternatively, intermolecular energy can be represented by van der Waals potential combined with DMA for electrostatic interactions.22g In that case, atomic multipoles can be derived directly from calculated molecular wave functions. A few examples of the successful application of the current CSP methods to predict structures of simple flexible molecules (up to three rotatable bonds) were reported in the literature (Figure 2). That includes carbamazepine,36 glycol and glycerol,37 aspirin,20 phenobarbital,38 piracetam,39 6-amino-2phenylsulfonylimino-1,2-dihydropyridine (compound VI second blind test),14b,21c,d alanine, valine, leucine and isoleucine amino acids,40 and 1,2-dichloro-4,5-dinitrobenzene (compound XVII fifth blind test).14e,41 In all these cases the experimental structures either corresponded to the global minimum in lattice energy or were among the lowest energy structures. A higher degree of conformational degrees of freedom, which is common for pharmaceutical compounds, would require a rational selection of preferred conformations for CSP study. An extreme example of a successful application to a flexible molecule crystal structure was recently reported for a benzyl(4(4-methyl-5-(p-tolylsulfonyl)-1,3-thiazol-2-yl)phenyl)carbamate molecule with 8 conformational degrees of freedom (compound XX fifth blind test, Figure 2m).14e,41 The two successful strategies (Day, Cruz-Cabeza and Galek; Price, Kazantsev, Karamertzanis, Adjiman and Pantelides) allowed the search space to be reduced to a more manageable level by taking advantage of different CSD database guided selections of starting conformations for the CSP studies.41 Both approaches predicted the observed crystal structure of compound XX as a global minimum in lattice energy. Another challenging application was reported for axitinib (Figure 2n),42 an oncology drug developed at Pfizer that inhibits vascular endothelial growth factor (VEGF).43 Axitinib is a particularly interesting compound with six rotatable bonds, which has been observed to form five anhydrous conformational polymorphs.44 To perform rational conformer selection, the OPLS_2005 force field 45 was modified to allow descriptions of close intramolecular contacts between a carbonyl oxygen and a sulfur (S···O) (the S···O distance in observed axitininib polymorphs is in the range 2.8−3.4 Å). The modified force field parameters allowed for the generation of preferred conformations that were accurately repacked into three out of four single conformation polymorphs, including the most stable observed form of axitinib (form XLI).42 Each of the computationally predicted forms was ranked among the top few empirically observed solutions. A novel approach for rational selection of preferred conformations for the CSP study of flexible molecules was
energy contributions to entropies and relative free energies of polymorphs of pharmaceutically related molecules have been considered in multiple studies. Van Eijck demonstrated that the rankings and (free) energy differences for glycol and glycerol improve upon introducing the corrections due to lattice vibrations.28 It was found that entropy and zero-point energy give the largest contribution to free energy differences between hypothetical crystal structures, adding up to about 3 kJ/mol for the structures with the lowest energy. The ranking of predicted crystal structures of pyridine improved upon inclusion of vibrational contributions to free energy.29 The predicted variation in the vibrational energy terms between the different structures is up to about 2 kJ/mol. In an earlier study, Gavezzotti and Filippini calculated vibrational contributions to entropy differences between 204 pairs of known polymorphs of organic molecules at 300 K.30 It was found that in the majority cases lattice vibration contributions to the differences between polymorphic free energies do not exceed 2 kJ/mol though it may approach a value of 5 kJ/mol. Therefore, based on these limited studies, an error in polymorph ranking based on the lattice energy landscape at 0 K, on average, may achieve up to 3 kJ/mol due to neglect of zero-point energy and entropic contributions at ambient temperature. Kinetics Effect. Another intrinsic limitation of polymorphism prediction by typical CSP techniques is related to the influence of kinetics on the observed crystal forms. It is reasonable to speculate that kinetics effects may also account for a much lower number of observed polymorphs relative to the predicted ones. Unfortunately, the current level of theoretical development and computational capabilities do not allow a detailed simulation of the crystallization process of pharmaceutical compounds. Therefore, simpler kinetic considerations were suggested and adopted in selected CSP studies, one of them being the growth rate of competing crystal forms based on a simple growth morphology theory.31 Initial studies allowed distinction between real and unobserved crystal structures of paracetamol, pyridine, and parabanic acid.7,29,32 Unfortunately, such important effects as solvent, degree of supersaturation, and temperature are ignored in a simple growth morphology theory, potentially limiting the applicability of this approach. An alternative approach, effectively taking into account kinetic effects, is based on a combination of an energy-based CSP prediction with a knowledge based approach based on Cambridge Structure Database (CSD)33 analysis of observed crystal structures.34 In fact, over half a million observed crystal structures in CSD effectively reflect valuable data set on both kinetic and thermodynamic preferences for crystallization. It is assumed (though is not proven) that the majority of crystal structures in the database correspond to stable forms. Molecular Flexibility and Conformational Polymorphism. Pharmaceutical compounds are typically quite flexible and, therefore, are subject to a conformational polymorphism. However, conformational polymorphism may be a challenging task for CSP study due to the following reasons. First of all, special attention should be given to an accurate and balanced treatment of intra- (conformational) and intermolecular energy contributions to the lattice energies of conformational polymorphs. In addition, the current CSP methods heavily rely on a correct selection of the starting molecular conformations that are typically held rigid at least during initial crystal packing generation. Since a drug-like molecule may display up to 10 rotatable bonds,35 a CSP study based on the 475
dx.doi.org/10.1021/op300274s | Org. Process Res. Dev. 2013, 17, 472−485
Organic Process Research & Development
Review
proposed by Abramov et al.46 Since crystallization in the pharmaceutical industry is not performed from the gas phase but rather from different solvents, COSMO-RS theory47 application was proposed for a conformational distribution analysis to select a set of preferred conformations for a virtual polymorph screening. In support of the approach, a correlation between trends of predicted (and, in some cases, also observed) changes of conformational distribution in different solvents and crystallization results was demonstrated for four cases of conformational polymorphs: n-phenylhydroxamic acid, ibuprofen, famotidine, and ritonavir. A small diverse set of solvents was proposed to be considered for selecting molecular conformations for virtual polymorph screening through CSP: polar protic water (both H-bond donor and acceptor capabilities) and diethyl amine (H-bond donor capabilities); polar aprotic acetone (H-bond acceptor capabilities); nonpolar hexane; and self-media to mimic solid amorphous. A combined set of conformations, which displayed the highest populations in any of the above solvents, was recommended for the virtual conformational polymorph screening. This conformational selection approach has been successfully applied for several compounds, one of which is the potent and selective dual ALK and c-MET inhibitor crizotinib (5 rotatable bonds, Figure 2o) developed for the treatment of ALK-positive nonsmall cell lung cancer. The CSP analysis results for crizotinib are considered below in section 5. A similar approach may be applied to assess likelihood of conformational polymorphism.46 Number of Molecules in Asymmetric Unit. The prevalence of pharmaceutical multicomponent crystals, which by definition display more than one molecule in the asymmetric unit, was demonstrated in a recent analysis of a series of 245 polymorph screens.2 In fact, salt formation and cocrystallization are being increasingly used to improve the physicochemical properties of solid pharmaceuticals, including solubility, dissolution rate, bioavailability, and chemical and physical stability.2,48 Though there are multiple reports with a successful prediction of crystal structures of multicomponent crystals (cocrystal,11b,49 salt,11b,50 solvate,51 and hydrate52), such systems are computationally much more expensive and are subject to the significantly increased impact of all the limitations described above. As a result, the general effectiveness and reliability of the current CSP methods tend to quickly decrease with an increase in system complexity. Therefore, despite the prevalence and importance of pharmaceutical multicomponent crystals they represent a major challenge in the exploration of polymorphism by the CSP approach. Summary. Computational approaches to crystal structure and polymorph prediction have demonstrated significant improvements and method developments in recent years. However, due to the outlined limitations, it is still hard to estimate the reliability of the predictions for a drug-like molecule by the lattice energy minimization alone. Perhaps only the crystal energy landscape, in which the observed form is separated from all other structures (Figure 1a) or from a global minimum in energy (Figure 1c) by a substantial energy gap (≥4−5 kJ/mol), may provide a reliable insight into the risk assessment of the solid form selection. Therefore, for general computational support of solid form selection in the pharmaceutical industry, the application of CSP technology in combination with the other methods described below is greatly preferable.
3. HYDROGEN BONDING PROPENSITY CONSIDERATION An alternative simplified approach to the classification of polymorph stability is based on selected dominant interaction analyses. As an example, a good polymorph stability ranking of the systems, which are dominated by nonspecific van der Waals interactions, can be provided by the crystal density rule.53 However, pharmaceutical small molecules typically contain a certain number of donors and acceptors and therefore display H-bonding interactions in the solid forms.35 Since H-bonds are the strongest and most specific (directional) interactions, they typically play a dominant role in the crystallization (nucleation and growth) and stability of pharmaceutical solids.54 As a result, systems with H-bonding interaction are known to violate the density rule. This inconsistency is addressed in two approaches, developed for the ranking of organic crystal stability based on H-bonding propensity analysis. Knowledge-Based Model. A probabilistic approach to analyses of organic crystal stability was recently developed55 based on statistical analysis of hydrogen bonds in the CSD. This knowledge-based method analyzes potential H-bonding functional groups of a molecule to make predictions on statistically likely pairwise H-bonding interactions in the stable crystalline form. The polymorphic stability classification by the model is based on the rationale of best donor/best acceptor pairing and the corresponding hierarchy of observed and absent H-bonds in a particular case.56 Based on this knowledge-based approach, a structure that adopts the highest propensity Hbonds displays a low likelihood of the existence of a more stable form thus indicating that the polymorph screening, which had not discovered any other stable form, was likely to be complete. The presence of only low(er) propensity H-bonds, in contrast, indicates a high risk of the existence of a more stable form and that an additional study is likely to be required. The current methodology is implemented in Materials Mercury57 and involves five stages: selection of appropriate H-bonding functional groups in the molecule(s) in an asymmetric cell; training set selection based on CSD (and inhouse database, if available) data sampling; fitting of logit Hbonding propensity (LHP) model, using qualitative (existence or absence of specific H-bonding observations in the training set) and optional quantitative (such as competition, donor and acceptor aromaticities, and steric densities) descriptors;55 model validation on a test set; and model application to the crystal of interest. Applications of the knowledge-based approach to four systems were reported in the literature: ritonavir (Figure 3a),56 N2-(indol-3-ylacetyl)-L-asparagin,55 and lamotrigine.58 The H-bonding propensity analysis for the latter case was used to assess a potential cocrystal formation and is out of the scope of the present review. Propensity analysis of ritonavir showed that, in contrast to initially marketed form I, all four Hbonds in the discovered later stable form II display a propensity, π, higher than 0.32 (Table 1). The propensity of the strongest H-bond in form II, which is absent in form I, was found between the amide donor and hydroxyl acceptor (π = 0.551, Figure 3a, Table 1). Therefore, it was demonstrated that if the knowledge-based methodology had been available, then the H-bonding analysis could have helped, recognizing the high likelihood of a missing stable form (II) at an early stage of drug development. 476
dx.doi.org/10.1021/op300274s | Org. Process Res. Dev. 2013, 17, 472−485
Organic Process Research & Development
Review
pyrazole amine and amide oxygen (Figure 3b). However, this H-bond was observed in all known axitinib crystal structures, including stable form XLI, and thus does not allow a stability ranking of the five observed anhydrous forms. This example outlines a limitation of H-bonding propensity application to the polymorph stability ranking, which is addressed in more detail below. The rest of the axitinib donors and acceptors should display a noticeably weaker hydrogen-bonding propensity according to these calculations and are readily assessable for the formation of multiple solvates.44 Theoretical Approach. The results of knowledge-based calculations may depend on the availability of observations for certain H-bonding functional groups in the CSD (and inhouse) database. That is why, in order to obtain good statistics for model fitting, some functional groups are often presented in a generalized form, resulting in the loss of a specific nature of the H-bonding centers in the molecule of interest. Therefore, an independent and relatively simple theoretical approach was needed, which would treat all H-bonding features explicitly to allow a good ranking of H-bonding propensities. It should also be recognized that such an approach will inevitably introduce its own limitations related to the applied level of theory. In order to address this issue, Abramov60 performed a study of the correlation of two independent schemes of atomic charge calculations (electrostatic potential (ESP) fitted charges, Qesp, and polarization σ H-bonding charges) with Abraham’s experimental H-bonding scale61 and with some available results from theoretical calculations, based on complexation reaction consideration.62 Qesp are partial atomic point charges, which are appropriate to better reproduce the molecular electrostatic potential at the molecular surface.63 These charges are quite popular. For example, they are used in the AMBER force field64 and may be applied to describe the electrostatic energy contribution in Dreiding65 and Williams 9923 empirical force fields. The σ H-bonding charges (also defined as a first atomic H-bonding moment, Mx1,HB)47 are based on the acceptor or donor atomic polarization charge density in the infinite dielectric above or below, respectively, a threshold value. These are surface charges, which directly reflect both polarization and surface accessibility effects. It was demonstrated that a nice ranking of H-bonding interaction energies may be performed by a product of the σ (rather than Qesp) Hbonding charges of the interacting donor and acceptor atoms.60 As an example of that study, correlations of theoretical H-bond interaction energies between substituted phenols and neutral or ionic probes with phenolic donor, (O)H, and acceptor, O, charges are presented in Figure 4. Qesp atomic charges provided a relatively poor ranking of all considered intermolecular interactions (Figure 4a,b). The best correlations were found between the σ charges of the acceptor and donor and the interaction energies of substituted phenols with the ionic probes (Figure 4c,d), for which H-bonding contributions are expected to be the highest. Klamt et al. recently reported a systematic study, which demonstrates a very good correlation of 2465 DFT/COSMO H-bond energies with a product of polarization σ charge densities of donor and acceptor atoms.66 The fitted relationship yielded an RMSD of only 0.36 kcal/mol. This result provides an excellent justification of the preferred applicability of COSMO-RS theory for H-bonding propensity ranking. The resulting model equation has produced scaling factors for σcharge densities of different donors and acceptors, which can
Figure 3. H-bonding functional groups considered in LHP model of ritonavir (a)56 and axitinib (b).59
Table 1. LHP propensity, π, predictions for potential donor−acceptor combinations in ritonavir (as labeled in Figure 3a), and observed H-bonds in two polymorphic forms; data taken from Galek, Allen, Fábián, and Feeder56 donor
acceptor
π
±a
form I
form II
amide amide carbamate hydroxyl amide amide carbamate hydroxyl carbamate hydroxyl carbamate hydroxyl ureido ureido ureido ureido amide amide carbamate hydroxyl carbamate hydroxyl
carbamate hydroxyl carbamate carbamate amide ureido hydroxyl hydroxyl amide amide ureido ureido carbamate hydroxyl amide ureido thiazoyl A thiazoyl B thiazoyl A thiazoyl A thiazoyl B thiazoyl B
0.618 0.551 0.538 0.537 0.501 0.499 0.47 0.469 0.42 0.419 0.418 0.417 0.319 0.263 0.225 0.224 0.152 0.142 0.115 0.114 0.107 0.106
0.094 0.052 0.09 0.09 0.055 0.072 0.078 0.037 0.083 0.045 0.088 0.058 0.086 0.041 0.04 0.044 0.054 0.05 0.044 0.039 0.041 0.036
− − X − X − − − − − − − − − − X − − − X − −
− X − − − − − − X − − X X − − − − − − − − −
a The error bars of the propensities predicted within confidence interval of 95%.
There are multiple applications of knowledge-based approaches to solid form risk assessment in the pharmaceutical industry, among which two examples will be discussed in the current review: axitinib59 and crizotinib. Analysis of crizotinib is presented in section 5. In the case of axitinib, the strongest donor−acceptor pair was predicted to occur between the 477
dx.doi.org/10.1021/op300274s | Org. Process Res. Dev. 2013, 17, 472−485
Organic Process Research & Development
Review
Figure 4. Correlations of theoretical H-bonding interaction energies62 between substituted phenols and neutral or ionic probes with the charges of the H-bonding centers on phenolic hydroxyl (O acceptor and (O)H donor). Energies are presented relative to the unsubstituted system. ESP charges are used to describe interactions of substitutes phenols with protonated methylamine and acetate anion probes in (a) and (b), respectively. H-bonding interactions of substituted phenols with methanol (blue line) and protonated methylamine (red line) donor probes are ranked by σ acceptor charges in (c). H-bonding interactions of substituted phenols with formaldehyde (blue line) and acetate anion (red line) acceptor probes are described by σ donor charges in (d). ESP charges were calculated by DMol3 software81 at the PBE/DNP level of theory.82b σ-Donor and acceptor charges were calculated by COSMOtherm76 software using polarization molecular surfaces predicted at the PBE/DNP/COSMO level of theory.82
also be applied to the corresponding σ atomic charges for a better ranking of the H-bonding propensities. An application of H-bonding propensity analysis based on σ charges to support the solid form selection of axitinib was reported by Campeta et al.44a The largest positive and negative σ H-bonding charge values correspond to the strongest acceptor and donor features, respectively (Figure 5a). According to these results, the strongest donor−acceptor pair for this molecule should be the pyrazole amine−amide oxygen H-bond. This result is in excellent agreement with the predictions based on the knowledge-based model described above. σ H-bonding charges of ritonavir suggest that the strongest H-bonding pairing should take place between the ureido acceptor and hydroxyl donor, followed by the amide acceptor and carbamate donor (Figure 5b). These H-bonds are present in stable form II, while are absent in metastable form I. Therefore, the σ H-bonding charge analysis properly indicates the high risk of form I selection as a commercial solid form of
ritonavir. It should be noted that H-bonding propensity ranking based on σ charges disagrees with the results of the reported knowledge-based model,56 though both approaches indicate a high likelihood of a new stable form. 5-Methyl-2-[(2-nitrophenyl)amino]-3-thiophenecarbonitrile (ROY) is a famous nonpharmaceutical example of multiple polymorphs, existing at ambient temperature.67 The σ Hbonding charge analysis (Figure 5c) demonstrates that cyano nitrogen is the strongest acceptor in the molecule, and intermolecular H-bonding shows a preference for occurring in the stable form of ROY. That is supported by experimental observation that the most stable yellow polymorph (Y) is indeed stabilized by such a H-bond in contrast to metastable YN, ON, OP, and R forms. Application of the σ H-bonding charge analysis to support the commercial solid form selection of crizotinib is presented in section 5. Consideration of Limitations of Hydrogen Bonding Propensity Approach. There are several important limi478
dx.doi.org/10.1021/op300274s | Org. Process Res. Dev. 2013, 17, 472−485
Organic Process Research & Development
Review
Figure 5. Examples of theoretical H-bonding propensity analysis for axitinib (a), ritonavir (b), and ROY (c). The σ H-bonding charges were calculated at the BP/TZVP/COSMO level of theory. Acceptor and donor charges are highlighted in red and blue colors, respectively. Charges for the strongest acceptor and donor are displayed in bold.
fifth blind test was performed based on both the lowest lattice energy and LHP model considerations.41 Alternatively, the results of the initial crystal packing generation, which do not display the preferred strongest H-bonding interactions, could be filtered out from the following CSP analysis. In that case, the final polymorph ranking will fully satisfy H-bonding propensity expectations and will be based on the lattice energy minimization results. In an attempt to overcome the above limitation, Abramov suggested taking an extra step beyond H-bonding propensity analysis and focused on the estimation of the relative strength of crystallographically observed intermolecular H-bonds in different polymorphic forms of drug-like molecules.72 It was demonstrated that prediction of the most stable form may be performed based on the relative strength of the strongest Hbonding interaction, as measured by QTAIM (Quantum Theory of Atoms in Molecules) 73 properties at the corresponding bond critical point. For example, a linear relationship between ΔHfus data and the potential energy
tations in stability ranking based on H-bonding propensity analysis in addition to the fact that H-bonding is not the only contribution to the stability of drug-like crystals. Similar to the CSP approach, the models based on the propensity of the Hbonding interaction generally cannot account for the enantiotropic relationship between the polymorphic forms at the temperature of interest (typically ambient temperature). In addition, an inherent limitation of the knowledge-based method is the inability to estimate the stability ranking of the polymorphs displaying similar H-bonding patterns in the solid state. There are several known examples of the drug-like polymorphic systems, which preserve at least the strongest Hbonding interactions for two or more forms. Among them are all five axitinib anhydrous forms,44a rotigotine forms I and II,68 paracetamol forms I and II,69 nimodipine forms I and II,70 and four carbamazepine anhydrous forms.71 One of the ways to address this issue is to combine a CSP study and H-bonding propensity analysis. For example, the final selection of three predicted structures of compound XX in the 479
dx.doi.org/10.1021/op300274s | Org. Process Res. Dev. 2013, 17, 472−485
Organic Process Research & Development
Review
density, Vb, at the bond critical point between the pyrazole amine and amide oxygen was found for four axitinib polymorphs (Figure 6). The actual rank order of the
amine and amide oxygen. This observation demonstrates the dominant contribution of the strongest H-bonds to relative lattice energies of the related polymorphs. The proposed QTAIM approach may help, predicting the relative stabilities of drug polymorphs from experimentally observed or computationally predicted crystal structures. Alternatively, qualitative consideration of H-bonding networks may be applied to allow an estimation of the stability ranking of the polymorphs displaying similar H-bonding patterns in the solid state. For example, for all axitinib anhydrous forms except for form XLI, the strongest H-bond between the pyrazole amine and amide oxygen leads to the formation of molecular dimers, which are connected to each other by additional weaker H-bonds (Figure 7a−d).44a The strongest H-bonds form an extended network throughout the crystal only in form XLI (Figure 7e). It was recognized by Campeta et al.44a that the presence of this strong H-bond network is favorable for the thermodynamic and mechanical stability of the crystal. For example, the destruction of the long order packing through melting of all forms, except for form XLI, can be achieved by imposing a certain amount of energy, sufficient to break only the weak H-bonds connecting the molecular dimers. Greater energy is required to break the chains of the strongest H-bonds in form XLI.
Figure 6. A correlation between Vb descriptor as a measure of relative strength of H-bond between pyrazole amine and amide oxygen and ΔHfus properties of axitinib polymorphs.72 Reproduced from ref 72 with permission.
4. SOLVENT SELECTION FOR STABLE FORM CRYSTALLIZATION A common method of most stable form discovery is slurry crystallization. In this experimental approach, a metastable form
experimental stabilities of the anhydrous forms at ambient temperature, XLI > XXV ≥ VI > IV > I,44a is reasonably well described by the strength of the H-bond between the pyrazole
Figure 7. Crystal packing patterns of axitinib anhydrous forms: (a) form I, (b) form IV, (c) form VI, (d) form XXV, and (e) form XLI.44a The strongest H-bonding between the pyrazole amine and amide oxygen takes place in all the forms; however only in the stable form XLI (e) it forms chains throughout the crystal. Reproduced from ref 44a with permission. 480
dx.doi.org/10.1021/op300274s | Org. Process Res. Dev. 2013, 17, 472−485
Organic Process Research & Development
Review
Table 2. Results of polymorph screen of ritonavir74 and predicted conformer populations at ambient temperature;a data taken from Abramov, Zell, and Krzyzaniak46
a
solvent
solid form (2 days)74
solid form (2 weeks)74
conformer I population, %
conformer II population, %
water hexane MTBE 1,2-xylene toluene nitromethane ethyl acetate acetonitrile 2-propanol 2-butanone acetone 1,2-dimethoxyethane ethanol
I I I I I II II II II II II II II
II I,II mixture I I I II II II II II II II II
21 85 39 70 68 31 35 26 32 34 27 28 30
79 15 61 30 32 69 65 74 68 66 73 72 70
Predictions by COSMOtherm software at the PBE/DNP/COSMO level of theory.46
was demonstrated that, except for MTBE, the results of the 2week slurry crystallization follow the trend of the preferred conformer population in the corresponding solvent. Therefore, it was proposed that the predicted trend of conformer population in different solvents may allow tuning (promotion or inhibition) of the corresponding conformational polymorph crystallization. The final selection of solvents should also satisfy the criteria of relatively high (predicted or measured) solubility of the API. This approach was applied to a targeted polymorph screening of axitinib.44a As discussed above, axitinib conformations which may form an extended network of the strongest H-bonds, like that in form XLI, would be expected to be significantly more stable than dimer-based structures observed in all other forms. It was shown that, according to the results of a full conformational search, other axitinib conformations (‘‘form XLI-sister’’ conformers) may exist which were not observed in the early polymorph screening, but could also potentially produce forms with an extended network of the strongest H-bonds. It was proposed that a significant increase of the experimental temperature may result in breaking up axitinib dimers, increasing the populations of unobserved form XLI-sister conformers. For validation of this rationale, the temperature dependence of the conformational populations of axitinib in a number of high boiling solvents was calculated, using COSMO-RS theory. The calculations demonstrated an increase of the absolute and relative populations of the form XLI-sister forms with temperature in all solvents. However, the population of the form XLI conformer stayed noticeably higher than the population of any of its sister forms, favoring formation of the form XLI polymorph. Based on the calculations, several high boiling solvents were selected for an experimental followup. No new forms of axitinib were observed from this study. Therefore, the results provided confidence that the current development form was the most stable polymorph, with a low likelihood for the existence of a more stable anhydrous form. Another important task in the pharmaceutical industry is the solvent selection for stable form screening, which does not tend to form solid solvates (sometimes called pseudopolymorphs) with the API. Hydrate formation can be considered as a particular case of this broader issue. Solvate formation is typically undesirable due to the following reasons. Solvates
is suspended in a saturated solution and recrystallized into the stable form via solvent mediated transformation. It was demonstrated that the transformation rates are faster in solvents which provide higher compound solubility; a solubility of at least 8 mM (but
