Pharmaceutical co-crystals, comprising an active pharmaceutical ingredient (API) and another pharmaceutically acceptable molecule, can be used to address intellectual and physical property issues in pharmaceutical development without changing the chemical composition of the API. The search for suitable co-crystals is hindered by the lack of computational tools to investigate whether the formation of the heteromeric is thermodynamically favoured over the homomeric crystals, and predict their structures.
This paper presents a novel methodology for the prediction of co-crystal structures, based only on the atomic connectivity of the component molecules and assumed stiochiometry. The method uses global crystal energy minimisation, for both the single- and two-component crystals. This allows a comparison of the relative stability of the co-crystal with respect to the energies of the two pure-component crystals, and can therefore predict when co-crystal formation is thermodynamically favoured.
Crystal structure prediction of co-crystals is a challenging problem because of the presence of two, often flexible, crystallographically independent molecules in the asymmetric unit. Hence, a search that is sufficiently broad to reliably locate all energetically competitive minima on the crystal energy surface can typically be performed only with empirical, computationally inexpensive models that may lack sufficient realism for accurate prediction of relative energies. Therefore, the results of such a search need to be refined using more accurate models for both the intra- and inter-molecular contributions to the lattice energy.
In this paper, we propose a 3-stage algorithm:
· Stage I: Global Search for low-energy minima.
This is carried out using an approach1 based on low-discrepancy sequences. The intramolecular energy is approximated using Hermite interpolants2 on a grid of isolated-molecule energies, and its gradients with respect to all torsion angles that can be affected by intermolecular forces. The grid is pre-calculated using a number of quantum mechanical energy minimisations. The intermolecular interactions are modelled using a computationally inexpensive isotropic model, comprising an empirical repulsion-dispersion potential3;4 and a set of conformationally invariant atomic charges.
· Stage II: Re-ranking of low-energy crystal structures based on improved electrostatic interactions.
The low-energy minima computed at Stage I are used as initial points for local energy minimisations. In this case, intramolecular energy and intermolecular repulsion-dispersion contributions are computed as at Stage I, but the intermolecular electrostatic interactions are modelled5 by a set of distributed multipole moments on each atom up to hexadecapole.6 The variation of the charge density with conformation is restricted to the analytic rotation of the atomic multipole moments with the local environment of each atom.
· Stage III: Re-ranking of low-energy crystal structures based on conformation-dependent electrostatic interactions and accurate intramolecular energies.
The low-energy minima computed at Stage II are used as initial points for local energy minimisations involving a nested 2-level optimisation approach.7 The outer level manipulates the molecular conformation (i.e. the flexible torsion angles) while the inner one determines the optimal crystal structure for a molecule of fixed conformation (as determined by the current outer iteration). The intra-molecular energy and the atomic multipole moments are determined at each outer iteration using “on-the-fly” quantum mechanical calculations. Given the cost of this calculation, Stage III is applied only to very few (lowest-energy) crystal structures determined by Stage II.
The above algorithm is applied to the prediction of the structure of the co-crystals of p-aminobenzoic acid with 2,2'-bipyridine. The application of the algorithm to pure p-aminobenzoic acid resulted in its low-temperature polymorph8 being ranked 1st i.e. being identified as the one of lowest energy. The experimentally known crystal structure9 of 2,2'-bipyridine was ranked 2nd. The known p-aminobenzoic acid/2,2'-bipyridine crystal structure10 is also ranked 1st. Since its energy is almost equal to the sum of the global minimum crystal energies of the component molecules, the predicted thermodynamic driving force for forming the cocrystal is small. We discuss the limitations of this methodology, including the possibility that p-aminobenzoic - 2,2'-bipyridine co-crystal formation may be kinetically rather than thermodynamically driven.
Keywords: co-crystals, crystal structure prediction, polymorphism, global optimisation
Reference List
(1) Karamertzanis, P. G.; Pantelides, C. C. J.Comput.Chem. 2005, 26, 304-324.
(2) Karamertzanis, P. G.; Pantelides, C. C. Mol.Phys. 2007, 105, 273-291.
(3) Williams, D. E. J.Comput.Chem. 2001, 22, 1154-1166.
(4) Coombes, D. S.; Price, S. L.; Willock, D. J.; Leslie, M. J.Phys.Chem. 1996, 100, 7352-7360.
(5) Willock, D. J.; Price, S. L.; Leslie, M.; Catlow, C. R. A. J.Comput.Chem. 1995, 16, 628-647.
(6) Stone, A. J. J.Chem.Theory Comp. 2005, 1, 1128-1132.
(7) Karamertzanis, P. G.; Price, S. L. J.Comp.Theor.Comp. 2006, 2, 1184-1199.
(8) Lai, T. F.; Marsh, R. E. Acta Crystallogr. 1967, 22, 885-&.
(9) Kuhn, F. E.; Groarke, M.; Bencze, E.; Herdtweck, E.; Prazeres, A.; Santos, A. M.; Calhorda, M. J.; Romao, C. C.; Goncalves, I. S.; Lopes, A. D.; Pillinger, M. Chemistry-A European Journal 2002, 8, 2370-2383.
(10) Smith, G.; Lynch, D. E.; Byriel, K. A.; Kennard, C. H. L. J.Chem.Crystallogr. 1997, 27, 307-317.