Computational polymorph screening reveals late-appearing and poorly-soluble form of rotigotine

The active pharmaceutical ingredient rotigotine—a dopamine agonist for the treatment of Parkinson’s and restless leg diseases—was known to exist in only one polymorphic form since 1985. In 2008, the appearance of a thermodynamically more stable and significantly less soluble polymorph led to a massive batch recall followed by economic and public health implications. Here, we carry out state-of-the-art computational crystal structure prediction, revealing the late-appearing polymorph without using any prior information. In addition, we predict a third crystalline form of rotigotine having thermodynamic stability between forms I and II. We provide quantitative description of the relative stability and solubility of the rotigotine polymorphs. Our study offers new insights into a challenging polymorphic system and highlights the robustness of contemporary computational crystal structure prediction during pharmaceutical development. It is important to screen potentially polymorphic systems during pharmaceutical development. Here, the authors demonstrate computational crystal structure prediction analysis which is capable of identifying the late appearing low-solubility polymorph of rotigotine without prior information.

D eveloping crystalline forms with desired physico-chemical properties 1,2 is of significant importance for the pharmaceutical industry. Yet this family of crystalline solids is susceptible to the well-established phenomenon of polymorphism originating in their complex energy landscapes with many closelying local minima 3 . Most of these experimentally obtainable crystal-packing arrangements differ in their stability by less than 2-4 kJ mol −1 4 . Despite the advances in crystallization techniques, engineering molecular crystals during the drug development process is a non-trivial task. Often unexpected polymorphs emerge due to a complex interplay between thermodynamics and kinetics. The disappearance of the previously known stable form or sudden appearance 5 of a more stable form can have severe consequences. For instance, the new form may present properties that are not suitable for the intended purpose leading to significant economic and public health repercussions 5,6 .
One recent example of a late-appearing polymorph is rotigotine, the active substance of Neupro®-a transdermal patch used in the treatment of Parkinson's disease (PD) and restless legs syndrome (RLS) [7][8][9] . It is estimated that about 6.3 million people worldwide are living with PD while RLS has a prevalence of about 7-10% of the adult population worldwide. Rotigotine is a dopamine agonist, stimulating the brain so that patients can control their movement and have fewer signs and symptoms of Parkinson's disease, such as stiffness and slowness of movement 9 . The delivery system consists of a fully amorphous patch formulation made of a molecular dispersion of rotigotine within a polymer matrix. Until 2007, the only and long-known crystalline phase of rotigotine was form I. In 2008, however, the promising clinical use of Neupro® was interrupted after the appearance of snowflake-like crystals within the patch 10,11 . This quality defect led to a massive batch recall in March 2008 in Europe and an outof-stock situation in the United States. The unexpected crystallization happened due to the appearance of a novel solid phase with greatly enhanced thermodynamic stability. The new solid was characterized as a conformational polymorph of rotigotine and subsequently named form II 10 . The solubility of form II was found to be more than eight times lower compared with form Ia striking difference for polymorphism, where mostly solubility ratios below two are observed 3,12 .
The consequences for the patients were dramatic as, for instance, Neupro® patches remained unavailable in the United States from 2008 to 2012, until the patch was properly reformulated to a stable amorphous matrix. The current pharmaceutical composition is devoid of any crystalline defects and ensures therefore its therapeutic efficacy, despite the existence of form II. The discovery of a second crystalline polymorph is especially astonishing for rotigotine as it is a commercial drug that was known since the mid-eighties and was well investigated for decades before the discovery. No indication of the presence of a second crystalline form was ever observed in the early years of formulation development despite a polymorphism screening done in the early 2000s. The antiviral compound ritonavir, which is in fact the most cited example of disastrous polymorphism in the pharmaceutical industry, shows a remarkable resemblance with rotigotine since also here the polymorphism is conformational in nature. Moreover, the new form has a much higher stability and hence much lower solubility with an instability of the amorphous formulation as the most striking effect 13 .
The likelihood of such a late-appearing polymorph can nowadays be assessed via computational molecular crystal structure predictions (CSP), which yield the thermodynamical stabilities for a variety of different crystal-packing arrangements 14 . Computational CSPs provide an unbiased thermodynamical landscape in contrast to experimental screenings. In several cases these computational predictions have led to the subsequent experimental discovery of a new polymorph [15][16][17] . However, these predictions face many challenges because of the high-dimensional nature of the crystallographic and conformational space of a given molecule. In addition, due to the very similar stability of the different polymorphs, relative free energy calculations with an accuracy of about 1 kJ mol −1 would be needed. Regularly, the Cambridge Crystallographic Data Centre (CCDC) organizes blind tests in order to determine the quality of a variety of different CSP methods. [18][19][20] . In recent years, several improvements have been made, especially for methods relying on density-functional theory (DFT), in order to cope with the increasing difficulty of these blind tests. We have recently shown that a combination of the most successful crystallographic space sampling from the latest blind test with a sophisticated firstprinciples energy-ranking approach yields excellent results across all target systems of that blind test 21 . In addition, this approach can be applied to large and complex systems of pharmaceutical relevance.
In the present study, we demonstrate computational CSP is capable of identifying the late-appearing polymorph of rotigotine without prior information from experiments. Furthermore, we show that our search led to the prediction of the new crystalline form or polymorph, which was not seen before, with thermodynamic stability lying between forms I and II. We outline a step-by-step guide to CSP and the subsequent reliable energy calculations, and provide detailed quantitative and qualitative insights of relative stability and solubility of rotigotine polymorphs. We highlight that the synergy between modern computational CSP and first-principles electronic-structure methods can be a viable approach to obtain invaluable insights into the drug development, and could have had a significant impact on the case of rotigotine.

Results
Crystal structure prediction of rotigotine. The polymorphic energy landscape for the 25 structures from step 3 of the applied CSP procedure (see Methods) is shown in Fig. 1b with a cutoff of 25 kJ mol −1 from the discovered global minimum. Each point corresponds to a local lattice energy minimum based on PBE + NP energies. The structures of the two experimental polymorphs have been published as single-crystal structures. The metastable structure with P4 3 space-group symmetry (form I) appeared as Private Communication to the Cambridge Structural Database in 2001 22 , whereas the stable structure with P2 1 2 1 2 1 space-group symmetry (form II) as part of a Patent in 2009 23 . In both singlecrystal structures, the thiophene rings randomly adopt one of the two alternative orientations, related by a 180°rotation, which has been modeled as disorder in the experimental structures. In the calculations, the number of degrees of freedom had been restricted to one molecule per asymmetric unit, so only fully ordered crystal structures, corresponding to all thiophene rings throughout the crystal adopting the same orientation, can be generated. Indeed, in the CSP study each experimental polymorph is found twice, once with all thiophene rings in one orientation and once with all thiophene rings in the alternative orientation. As a result, the two most stable predicted structures (blue circles) both correspond to form II, and the two structures describing the disordered form I (red squares) are found to be about 6 kJ mol −1 less stable than form II. We have included the structural overlays, which show perfect agreement for all four ordered structures, in the SI ( Supplementary Fig. 2). Calculations performed with larger supercells modeling explicitly disordered structures, in which different thiophene rings in the same crystal structure can adopt different orientations, confirmed that the orientations of the thiophene rings are randomly distributed over the two alternatives. Hereafter, we will always use the 50%:50% average of the two respective configurations to describe the disordered forms I and II.
New crystalline form of rotigotine. The CSP also predicted a new form with P3 2 space-group symmetry (green triangle), denoted as predicted form III, lying between form I and form II in the crystal landscape. Since the stability of form III falls between the two experimentally observed forms, this structure is a potential candidate for crystallization experiments. It is not obvious from the comparison of the crystal structure of predicted form III to the structures of forms I and II, how form III may be crystallized. A successful crystallization experiment would have to both inhibit the appearance of the substantially more stable form II and promote the appearance of the predicted form III. All other predicted structures are less stable than the experimentally observed forms. We also observed that the number of structures considered in steps 1 and 2 for rotigotine is much smaller than the typical values shown in Fig. 1a. This is a common feature for compounds that have only a few highly stable packings separated from the other structures by a substantial energy gap. In such cases, only a few structures fall into the target energy window, even in step 1.
At a time where form II had not yet been discovered, the energy landscape shown in Fig. 1b, would not only have flagged rotigontine as a compound with a missing, substantially more stable crystal form, but also have provided a clear indication for how to obtain that missing form. Form II is predicted to have significantly higher density than form I. By crystallization from solution under pressure (see ref. 16 ), the energy gap between forms I and II can thus be increased further. In such an experiment it may be possible to cross the metastable zone width of form II before the saturation concentration of form I is reached, resulting in the crystallization of form II.
Relative polymorphic stabilities of rotigotine. The discussed CSP approach provides an excellent sampling of conformational and crystallographic spaces, however, a reliable calculation of relative polymorphic stabilities requires using higher levels of theory than PBE + NP. Aiming at a robust quantum-mechanical description of free energies, we performed a final energy ranking of the 25 predicted structures utilizing the approach introduced in ref. 21 . This method employs a robust hierarchical first-principles approach based on DFT by incorporating key physical contributions including many-body van-der-Waals dispersion interactions (MBD), a sophisticated treatment of electron exchange and correlation via a hybrid functional, and vibrational free energies within the harmonic approximation. Furthermore, anharmonic effects to the free energy could be considered by accounting for thermal expansion via the quasi-harmonic approximation and the use of Morse oscillators as discussed in ref. 21 . We found that these effects can modify relative stabilities on average by about 1 kJ/mol 24 . Therefore, these effects should have an insignificant influence on the presented results since rotigotine's energy landscape is rather sparse. It was shown that the PBE0 functional 25 in combination with the MBD model 26,27 is able to describe relative stabilities for several molecular crystals within 1 kJ mol −1 compared with experimental observations 17,28-30 . Here we use this method for calculating static energies, whereas geometry optimizations and vibrational free energies were carried out at the PBE + MBD level 26,27,31 , leading to an improved description of vibrational free energies 24 . All calculations were performed using the all-electron code FHI-aims 32 Supplementary Fig. 4. It can be seen that the stability ordering of forms I, II, and predicted form III are identical, however, there are substantial differences in the relative energies. In contrast, we observe different stability orderings for higher-energy structures.
In order to assess the reliability and accuracy of the relative stabilities calculated by DFT + MBD, we experimentally measured enthalpy differences between forms I and II, and the solubility difference between the two forms. Figure 2a shows the cumulative enthalpy curves obtained from differential scanning calorimetry measurements (see the Methods section for more details). The measured enthalpy difference between form I and form II amounts to 7.5 kJ mol −1 as shown in the figure. The calculated relative stability of the three rotigotine polymorphs obtained by four methods of increasing sophistication is shown in Fig. 2b. It is apparent that the ranking based on PBE0 + MBD total energy is in better agreement with the experimental relative enthalpy than that of PBE + NP. However, these two methods do not include vibrational contributions. In order to directly compare our stability ranking with the relative enthalpy from the experiment, we add the harmonic vibrational internal energy U vib to the PBE0 + MBD lattice energies, yielding the PBE0 + MBD + U vib ranking. The U vib is defined as E ZPE þ R C V dT, where E ZPE is the zero-point energy while C V and T stand for the heat capacity and temperature, respectively. It is noted that a pV term is not included as its impact on the rankings at ambient conditions is negligible.
This results in a stability difference of 7.6 kJ mol −1 between form I and II, which is in excellent agreement with the experimental value of 7.5 kJ mol −1 . This result further suggests the significance of employing rigorous energy calculations after the initial CSP, as the relative enthalpy estimation by PBE + NP, used during CSP procedure, is 18% off from the experiment (see the CSP ranking in Fig. 2b).
Going from relative enthalpies to relative Helmholtz free energies at 300 K (see PBE0 + MBD+F vib ranking in Fig. 2b) does not affect our conclusions. The vibrational free energy (F vib ) contributions close the gap between form I and form II by more than 2 kJ mol −1 , while form III remains more stable than form I. In order to elucidate the origin of this free energy change, we calculated the relative PBE0+MBD+F vib energies by considering for F vib only contributions of the phonon density of states up to certain wave numbers (see Fig. 2c). The relative stabilities are normalized to the most stable structure (form II) and approach at high frequency the full PBE0+MBD+F vib results. One can observe that the relative stabilities are fluctuating in the lowfrequency region but already nearly converge to the final PBE0 +MBD+F vib values after considering vibrations up to about 160 cm −1 . This clearly shows that the relative free energies are dominated by low-frequency vibrations. The phonon modes in this frequency range consist mostly of intermolecular translations and rotations, for which an accurate treatment of van-der-Waals dispersion interactions provided by the MBD method is crucial. The low-frequency phonon densities of states are shown in Supplementary Fig. 5. All three polymorphs have distinct vibrational spectra in the THz regime. These differences in the low-frequency vibrations are utilized in THz spectroscopy to detect certain explosives and drugs, as well as to experimentally distinguish between polymorphs 33 .
Relative polymorphic solubility of rotigotine. A crucial property of an active pharmaceutical ingredient is its solubility, which can differ significantly between different polymorphs as it is linked to the free energy of the crystalline form. The measured solubility of form I is nearly 8.4 times larger than that of form II at room temperature in ethanol. The solubility ratio of two polymorphs in the same solvent can be calculated from their free energies as K I / K II = exp[−ΔF/RT], where ΔF is the free energy difference, and R and T are gas constant and temperature, respectively. This follows directly from the definition of the equilibrium constant, K = exp [−ΔF 0 /RT], and the observation that the free energy of the solution is independent of the polymorph. We outline the derivation of this equation in the Supplementary Discussion. Using the estimated −5.28 kJ mol −1 for the free energy difference at room temperature based on PBE0 + MBD + F vib (see Fig. 2b), we obtain a solubility ratio of 8.3 for form I/form II, which is in excellent agreement with the experimental measurements. The summary of our results along with experimental measurements are tabulated in Table 1.

Discussion
We have demonstrated that a contemporary computational CSP method reveals the late-appearing polymorph of rotigotine, and is  capable of correctly describing complex polymorphic energy landscapes of pharmaceutically relevant compounds. Furthermore, we predicted a potential new form of rotigotine with a relative stability between form I and II. The obtained result of computational CSP could, if they were available at the time of the drug development of rotigotine, have led to additional efforts to find the predicted stable form, thus avoiding the dramatic late appearance during production. In addition, we have shown that modern first-principles electronic-structure methods (the DFT + MBD framework) can be used to determine the relative enthalpies and relative solubilities of polymorphs in unprecedented agreement with experimental measurements (see Table 1). The computational description of the solubility is especially challenging since it does not only require the correct description of the enthalpy but also accurate vibrational entropies.
It has recently been shown in a metastudy over 41 commercial crystal structure predictions of pharmaceutical compounds that the thermodynamically stable form is expected to be missing in 15-45% of the cases 34 . The accurate results-within experimental precision-obtained for rotigotine with respect to both polymorphic enthalpy difference and relative solubility difference suggest that the DFT + MBD framework has the potential to confirm and narrow down this range, to pinpoint candidates for late-appearing polymorphs more clearly, and to assess the maximum relative solubility drop to be expected in case of a disappearing-appearing-polymorph event, thus providing crucial quantitative information for risk management when developing pharmaceuticals.

Methods
Crystal structure prediction within GRACE. The crystal structure prediction (CSP) for rotigotine was performed using GRACE 2.4 software package 35 without using any a priori knowledge from experiments. The procedure is discussed in detail in the Supplementary Information (see Supplementary Note 1) and consists of two parts (see Supplementary Fig. 1): (i) a tailor-made force field 36 is parameterized from scratch with the Force Field Factory module and (ii) the actual CSP is performed within the CSP Factory module 36 . The CSP Factory uses a three-step procedure for each Z′ value considered (Z′ is the number of molecules in the asymmetric unit) as shown in Fig. 1a. It is noted that all single bonds shown in Fig. 1a were treated as independent degrees of freedom in the crystal structure generation. First, crystal structures are generated with the tailor-made force field by a Monte Carlo parallel tempering algorithm (step 1). Promising candidates are subjected to a dispersion-corrected density-functional theory (DFT-D) lattice energy optimization. In this case the Perdew-Burke-Ernzerhof (PBE) 31 functional is used together with the Neumann-Perrin (NP) dispersion model 37 . These optimizations are performed with the Vienna Ab initio Simulation Package (VASP) 38-40 using loose convergence settings (step 2). Finally, the most promising candidates from step 2 undergo a second PBE+NP lattice energy optimization with tighter convergence settings (step 3). Statistical methods are used to control the completeness of the three-step procedure by looking at the standard deviations σ 1 and σ 2 obtained for the force-field energies of the first step and the coarse DFT-D energies of the second step, respectively.
All DFT-D calculations were carried out with the PBE functional 31 and the dispersion correction according to Neumann-Perrin 37 . DFT calculations use a plane wave cutoff energy of 520 eV and a k-point spacing of roughly 0.07/Å. All lattice energy minimizations of the final step have been converged to within at least 0.003 Å for atomic displacements, 0.00025 kcal mol −1 /atom for energy changes, 0.7 kcal mol −1 /Å for atomic forces and 1.0 kbar for cell stress. In the second step, lattice energies are converged to within at least 0.02 Å for atomic displacements, 0.001 kcal mol −1 /atom for energy changes, 7.0 kcal mol −1 /Å for atomic forces, and 15.0 kbar for cell stress.
Crystal structures were generated with one molecules per asymmetric unit (Z′ = 1). According to the statistics of the Cambridge Structural Database (CSD) 41 , 87.8% of the crystal structures of chiral molecules, as rotigotine is, contain one molecule per asymmetric unit. For molecules such as rotigotine, searches with two molecules per asymmetric unit (Z′ = 2) are technically feasible and often economically reasonable, but in this case the search range was deliberately limited to one molecule per asymmetric unit. For molecules of the size and flexibility of Rotigotine, it is because of the Z′ = 2 part of the calculation that CSP is perceived as an expensive computational technique. The CPU time cost of CSP grows exponentially with the number of degrees of freedom, and the cost increase for moving from Z′ = 1 to Z′ = 2 is roughly the same as for doubling the molecular size and flexibility at Z′ = 1. However, since most observed structures crystallize with only one molecule per asymmetric unit, many important results can already be obtained by searches with Z′ = 1 at a fairly acceptable cost, and the rotigotine study provides an impressive demonstration of this fact. The crystal structure generation was carried out in all chiral space groups. A value of Δ = 1.0 kcal mol −1 was chosen for the target energy window in which the completeness of the CSP procedure is statistically controlled. Typically, many of the considered structures fall outside this window and are reported in the final results.
The tailor-made force field was generated in 8 days, and actual CSP took 3 days. For post-CSP calculations, it approximately took 8 CPU hours to fully optimize a structure containing 188 atoms in the unit cell using PBE + MBD light setting, and extra 2 CPU hours for energy calculations using PBE + MBD tight setting on 6 nodes, where each node has 16 cores, altogether summing up to 2 days for the three main polymorphs, in addition to the CSP procedure (see Hardware details in the next section) 24 . Table 4 in Ref. 24 shows relative computational cost for one singlepoint energy and one force calculation for a particular molecular crystal containing 172 atoms in the unit cell, which allows an estimation of the total CPU time needed for the free energy based on various methods.
The CSP calculations were carried out on 48 Intel Xeon CPUs E5-2670 2.6 GHz 8 core (=384 cores) with InfiniBand network and 2 GB memory per core. The postCSP free energy calculations were carried out at DRACO extension cluster at the Max Planck Computing and Data Facility on Intel Haswell Xeon E5-2698 processors 2.3 GHz with InfiniBand FDR14 network and 128 GB memory per node.
Calculation of relative stabilities with the DFT + MBD framework. All DFT + MBD calculations were performed by using the all-electron code FHI-aims 32,42-47 . Full lattice and geometry optimizations were performed for all 25 polymorphs using PBE + MBD 26,27,31 with light species default settings for basis functions and integration grids. The number of k-points n in each direction was determined by the smallest integer for which n × a ≥ 28 Å, where a is the respective unit cell length. For the optimizations we used a convergence criterion of 0.001 eV/Å for force components.
For the stability ranking, we subsequently calculated the PBE0 + MBD energy as described in ref. 21 . First, the PBE + MBD energy is calculated with fully converged tight species default settings. Second, we add the difference between PBE0 + MBD and PBE + MBD calculated using light species default settings. This provides a good approximation of converged PBE0 + MBD energies at a much smaller computational cost 21 .
The vibrational free energies were calculated within the harmonic approximation by utilizing phonopy 48 and finite displacements of 0.01 Å. All necessary force calculations were performed with FHI-aims utilizing PBE + MBD with light species default settings. The vibrational free energy was evaluated within phonopy by using a fully converged q-point mesh.
Experimental determination of heat capacities and solubilities. Rotigotine Form I and Form II were provided by UCB Pharma SA (Anderlecht, Belgium). Ethanol was purchased from VWR (Darmstadt, Germany).
Solubility was determined by measurement of the rotigotine concentration of the supernatant of a crystal suspension by the following procedure. Suspensions of rotigotine in ethanol with a nominal concentration of about 1200 mg ml −1 in case of Form I and 200 mg ml −1 in case of Form II were prepared in 1 ml HPLC vials (Wicom, Heppenheim an der Bergstrae, Germany) and shaken overnight at 20°C (Thermomixer Eppendorf, Hamburg, Germany). The suspensions were then centrifuged at 14,000 rpm for 10 min at 20°C (Andreas Hettich GmbH, Tuttlingen, Germany). The supernatant solutions were subsequently sampled, diluted and subjected to HPLC analysis (Agilent 1100 VWD, Agilent, Waldbronn, Germany) for concentration determination.
Differential Scanning calorimetry (DSC) was performed on a Mettler DSC822 (Mettler Toledo, Greifensee, Switzerland) connected to a refrigerated bath at 10°C. Nitrogen purge flow was set at 50 ml min −1 . Analysis was carried out with a heating rate of 10°C min −1 in pierced aluminum crucibles, from 30 up to 140°C. The DSC curves are shown in Supplementary Fig. 6.
The relative enthalpy difference between the two forms was derived from the DSC signal by integration of the heat flow signal into enthalpy 49 . To this purpose the heat flow signal was first normalized with regard to the weight, then converted to heat capacity (C p ) by division by the heating rate. Knowing that the original DSC The relative enthalpy (ΔH), free energy (ΔF), and solubility between Form I and II of rotigotine as obtained by our computational approach compared with experimental measurements measurements were not calibrated for absolute C p determination and also knowing that we are only interested in relative differences between the two forms, the C p curves were adjusted (only by vertical shifting and slope correction) to coincide in the melt part of the curves above the melting point of Form II. In this temperature region both forms will end up in the same molten state, irrespective of the original crystalline form and will hence be in an identical state. Calculation of the enthalpy curves was done by setting both forms arbitrarily to 0 J g −1 at 140°C and consecutive integration of the C p curves from high to low temperature, thereby reducing the enthalpy content of both forms while decreasing in temperature. Following this procedure, the relative enthalpy difference between both forms could be determined as a function of temperature. See Supplementary Discussion for complimentary discussions.

Data availability
The data supporting the findings of this study are available within the article and its Supplementary Information. The structure files in the format of Crystallographic Information Files (CIF) and FHI-aims geometry are available in Supplementary Data 1.