Exploring the Formation of Multiple Layer Hydrates for a Complex

J. Phys. Chem. B 2009, 113, 5929–5937

5929

Exploring the Formation of Multiple Layer Hydrates for a Complex Pharmaceutical Compound Xin S. Zhao,† J. Ilja Siepmann,*,† Wei Xu,‡ Y.-H. Kiang,‡ Agam R. Sheth,‡ and Sami Karaborni‡ Departments of Chemistry and of Chemical Engineering and Material Science, UniVersity of Minnesota, 207 Pleasant Street SE, Minneapolis, Minnesota 55455, and Pharmaceutical Research & DeVelopment, Merck & Company Inc., P. O. Box 4, WP78-304, West Point, PennsylVania 19486 ReceiVed: September 13, 2008; ReVised Manuscript ReceiVed: January 18, 2009

The pharmaceutical compound A, 3-{2-oxo-3-[3-(5,6,7,8-tetrahydro[1,8]naphthyridin-2-yl)propyl]imidazolidin1-yl}-3(S)-(6-methoxypyridin-3-yl)propionic acid, is known to exist in five different crystalline forms that differ in the hydration state ranging from the anhydrous desolvate over hemihydrate, dihydrate, and tetrahydrate forms to the pentahydrate. The formation of the higher hydrates and the concomitant lattice expansion leads to undesirable tablet cracking at higher humidities. In this work, particle-based simulation techniques are used to explore the hydrate formation of compound A as a function of humidity. It is found that a simulation strategy employing Monte Carlo simulations in the isobaric-isothermal and Gibbs ensembles and transferable force fields, which are not parametrized against any experimental data for compound A, is able to yield satisfactory crystal structures for the anhydrate and pentahydrate and to predict the existence of all five hydrates. 1. Introduction In the pharmaceutical research and development process, a tremendous amount of effort and many bottlenecks separate a potential drug candidate with desirable biological activity from a product that can be given to patients. Formulation design is one of these bottlenecks because one needs to find an optimal formulation that is stable under various storage conditions for a reasonable period of time and has controlled bioavailability in the drug delivery process. Sometimes a drug candidate may possess physicochemical characteristics that make it extremely difficult to find an optimal or even acceptable drug formulation, but these difficulties are only fully understood after laborious experiments. With recent dramatic increases in the efficiency of computational algorithms, the accuracy of force fields, and computer speed, the research enterprise may reach a point where some of the initial screening of a drug candidate’s physicochemical properties can be carried out through in silico experiments. In addition, these computer simulations can provide molecular-level understanding of the formulations and their processing. Over the past decade, the prediction of the stable crystalline polymorph(s) for an organic molecule given only its chemical structure has received much attention.2-5 In contrast, the computational prediction of hydrate and solvate formation is in its infancy.1,2,6-10 On the other hand, molecular simulations have previously been used to understand the swelling of inorganic clays.11-14 The compound A, 3-{2-oxo-3-[3-(5,6,7,8-tetrahydro[1,8]naphthyridin-2-yl)propyl]imidazolidin-1-yl}-3(S)-(6-methoxypyridin-3-yl)propionic acid, intended for the treatment of osteoporosis, has recently been explored in the Merck Research Laboratories.15,16 This compound with the sum formula C23H29O4N5 exits in aqueous solution as a zwitterionic form (see Figure 1). In the drug development process, it was found that uncoated * Corresponding author. E-mail: [email protected]. † University of Minnesota. ‡ Merck & Company Inc.

Figure 1. Structure formula with atom index for compound A.

and coated tablets of compound A exhibit anisotropic volume expansion and develop macroscopic cracks when stored at elevated temperatures and humidities (see Figure 2). An experimental study of the moisture sorption isotherm revealed that compound A exists in five different hydration states including anhydrous, hemihydrate, dihydrate, tetrahydrate, and pentahydrate forms. The transition humidities range from about 15% relative humidity (RH) for the transition from anhydrous form to hemihydrate to about 80% RH for tetrahydrate to pentahydrate.15,16 The structure of the pentahydrate was determined from single-crystal X-ray diffraction (XRD), and those of the other hydrate forms were obtained from powder XRD. A comparison of the structures demonstrates that the compound A molecules form a layered framework and largely retain their conformation and relative position in these layers for all the hydration states, whereas the unit cell parameter perpendicular to the layers changes by ≈7 Å (≈20%) from the anhydrous form to the pentahydrate due to an expansion of the interlayer spacing. This expansion most likely causes the tablets to crack at high humidity. Another example of a crystalline pharmaceutical compound that exists in multiple hydration states is cromolyn sodium, which can absorb up to a stoichiometry of about nine

10.1021/jp808164t CCC: $40.75  2009 American Chemical Society Published on Web 04/07/2009

5930 J. Phys. Chem. B, Vol. 113, No. 17, 2009

Figure 2. Photo of an Avicel-based formulation of compound A stored at 313 K/75% RH for about 20 days.

water molecules and undergoes lattice expansion predominantly in one crystallogaphic direction.17,18 In this work, we use modern computational techniques to explore the hydrate structures and the adsorption isotherm for compound A. It is of particular interest to evaluate whether these computational techniques could be used to predict the existence of multiple hydrates. The remainder of this paper proceeds as follows. The simulation details are provided in section 2. Thereafter, in section 3, results for the structural properties and the moisture adsorption isotherm are presented, followed by some concluding remarks in section 4. 2. Simulation Details 2.1. Molecular Models. Compound A was represented by a united-atom model with parameters taken from the transferable potential for phase equilibria (TraPPE-UA) force field19-25 and the optimized potentials for liquid simulations (OPLS-UA) force field.26,27 These force fields utilize pseudoatoms located at carbon centers to model CHx groups, while treating all other atoms explicitly. The TraPPE-UA force field is parametrized to reproduce vapor-liquid coexistence curves over a range of temperatures, whereas the OPLS-UA force field is parametrized to reproduce the liquid density and heat of vaporization at or near the normal boiling point. Both force fields are fitted over a training set of low-molecular-weight compounds which contain the basic building blocks for larger, more complex compounds. The water molecules are represented by the popular TIP4P model.28 These force fields describe the intermolecular interactions through pairwise additive Lennard-Jones (LJ) and Coulombic interactions. The combination of TraPPE-UA and TIP4P models has been shown to accurately model octanol-water partitioning,29 retention in reversed-phase liquid chromatography,30 the solubility of CO2-philic polymers,31 and the phase behavior of benzene.32 Furthermore, the TraPPE-UA force field has also been applied to investigate the structure and adsorption isotherm of warfarin sodium 2-propanol solvate.9 For the present study, the intramolecular conformations of compound A and water were kept rigid during all simulations. The conformation of compound A was set to that determined for the pentahydrate from single-crystal XRD.15,16 Although compound A possesses multiple dihedral angles, it adopts very similar conformations in the five crystalline pseudopolymorphs investigated here as indicated by key dihedral angles and atom-atom distances (averaged over the distinguishable molecules in the unit cell). That is, the 1-12-13-14, 12-13-14-15,

Zhao et al. and 18-21-22-23 dihedral angles are 175°, -64°, and -66° for the pentahydrate, respectively, and the corresponding values for the anhydrous form are 175°, -57°, and -61°, respectively. The 7-23 and 7-30 distances are 13.30 and 10.76 Å for the pentahydrate, 13.31 and 10.61 Å for the tetrahydrate, 13.36 and 11.04 Å for the dihydrate, 13.10 and 10.70 Å for the hemihydrate, and 13.36 and 10.93 Å for the anhydrate. Hence, it appears to be a sensible approximation to treat compound A as conformationally rigid in the simulations. As indicated above, compound A exists as a zwitterion in aqueous solution, but the carboxylate group and the protonated nitrogen site 2 (see Figure 1) are in close proximity in the crystalline pentahydrate with two hydrogen bonds between the carboxylate group and the two nitrogen atoms (sites 2 and 4) on the tetrahydronaphthyridine ring. An exploratory simulation of the pentahydrate using a set of partial charges that would correspond to a complete transfer of the proton yielded rather disappointing structures, and one should expect only a partial transfer in this doubly hydrogen-bonded structure with similar charges on the two nitrogen-bound hydrogen atoms (sites 5 and 26). An additional test simulation where only the OPLS partial charge of a carboxylic acid proton, qH(COOH) ) +0.45e, was transferred to location 2 and the TraPPE partial charge for a hydrogen atom on a secondary amine was used for both hydrogens 5 and 26 yielded rather satisfactory results, and this set of partial charges was used for all subsequent simulations. Other partial charges were either taken directly from the TraPPEUA and OPLS-UA force fields or assigned to achieve neutral subunits of compound A. The complete set of LJ parameters and partial charges for the interaction sites of compound A is provided in the Supporting Information. A spherical potential truncation at 14 Å with analytical tail corrections33 was used for the LJ interactions. The LJ parameters for unlike interactions were computed using the standard Lorentz-Berthelot combining rules.34 The Coulombic interactions were computed using the Ewald summation method with tin foil boundary conditions35,36 and the following parameters: rcutoff ) 14 Å and κ ) 0.25. 2.2. Constant-Stress Monte Carlo Simulations of Anhydrate and Pentahydrate. As an initial step to validate the force field, constant-stress Monte Carlo (MC) simulations37-39 were used to investigate the crystal structures of the anhydrate and pentahydrate. The structure of the pentahydrate has been solved by single-crystal XRD,16 and it belongs to the C2 space group with two crystallographically independent compound A molecules in a unit cell, i.e., the numbers of compound A and water molecules per unit cell of pentahydrate are 8 and 40, respectively. Given the similarities of the powder patterns of the five pseudopolymorphs and the solved structure of the pentahydrate, the unit cell parameters of the other four forms could be deduced from their high-resolution X-ray powder patterns.16 The constant-stress MC simulations were carried out at T ) 298 K and p ) 1 atm. The simulation box for both the anhydrate and pentahydrate contained 80 compound A molecules (i.e., 2 × 5 × 2 and 1 × 5 × 2 unit cells for the orthorombic anhydrate and monoclinic pentahydrate forms) and an additional 400 water molecules for the pentahydrate. Phase space was sampled through rigid-body translations and rotations for compound A and water with the maximum displacement parameters set for each molecule type to achieve a 50% acceptance ratio. The volume and shape of the simulation box were allowed to fluctuate in contact with an external pressure bath. The simulation box was described by the H-matrix calculated from the six lattice parameters (a, b, c, R, β, and γ). Monte Carlo moves

Hydrate Formation in a Pharmaceutical Compound

J. Phys. Chem. B, Vol. 113, No. 17, 2009 5931

TABLE 1: Unit Cell Parameters of Anhydrous (A) and Pentahydrate (P) Forms obtained from Simulationa and Determined Experimentally16 form anhydrateb pentahydrate

method

a [Å]

b [Å]

c [Å]

R [deg]

β [deg]

γ [deg]

V [Å3]

sim exp sim exp

29.771 29.828 36.373 36.9615

7.601 7.4635 7.641 7.4582

20.431 20.6041 20.461 20.6914

90.01 90.00 90.01 90.00

89.92 90.00 102.82 99.461

95.01 90.00 90.32 90.00

46042 4587 55427 56264

a Subscripts denote the standard error of the mean in the last digit for the simulations and the uncertainty of the refinement for the experiment. b The experimental unit cell of the anhydrate is orthorombic and contains only four molecules of compound A. For better comparison with the monoclinic simulation cell and the pentahydrate form, the a cell parameter and volume are doubled.

were performed on the six nonzero elements of the H-matrix. While the experimental crystal structures of the anhydrate and pentahydrate forms are orthorombic and monoclinic, respectively, angular constraints were not imposed on the simulated systems. Both equilibration and production periods consisted of at least 105 MC cycles (one MC cycle consists of N randomly selected trial moves, where N is the total number of molecules in the system). 2.3. Moisture Sorption Isotherm Simulations. Gibbs ensemble Monte Carlo (GEMC) simulations40-43 with two different setups were used to investigate the hydrate states of compound A via moisture sorption isotherms. Initially, GEMC simulations with fixed unit cell parameters were carried out to verify that the simulations yield the correct hydration state for a given unit cell taken from the experimental single-crystal and powder XRD studies.15,16 Thereafter, a second set of GEMC simulations with variable unit cell volume and shape was initiated from the pentahydrate structure to explore whether the existence of the other hydrates could be predicted. A GEMC simulation utilizes two (or more) simulation boxes that represent different phases and are connected through special MC moves. GEMC simulations with a fixed unit cell for the first box containing the crystalline phase were performed on the hydrate structures that had been experimentally determined: hemihydrate, dihydrate, tetrahydrate, and pentahydrate forms. The second box contained a water vapor phase in contact with an external pressure bath that allows one to specify the appropriate RH. The GEMC simulation setup with fixed unit cell closely mimics a grand canonical MC simulation, but instead of specifying the chemical potential of an external reservoir, the former explicitly simulates the reservoir at the appropriate condition. This has the advantage that one does not need to know the equation of state of the reservoir molecules. The initial structure of the crystalline box contained 1 × 5 × 2 unit cells, i.e., 80 compound A molecules, using the unit cell parameters and molecular arrangement from the corresponding experimental structure15,16 but without any water molecules present. The vapor box contained initially 480 water molecules. The RH was increased stepwise in 10% intervals starting from RH ) 10%. That is, a given simulation was equilibrated for at least 105 MC cycles until there was no discernible drift in the average hydration number. This was followed by a production period of 0.5 × 105 MC cycles before increasing the RH to the next level. Both the sorption and desorption isotherms were computed. These simulations employed rigid-body translational and rotational moves, volume moves for the gas-phase box,37 and exchanges of water molecules using the configurational-bias MC technique.44-46 The latter move allows guest molecules to jump directly from the interior of the crystal to the vapor phase and vice versa without the need for an explicit pathway through which they could diffuse. The GEMC simulations with variable unit cell were started from the pentahydrate. The desorption isotherm was followed

by decreasing the RH stepwise from RH ) 100% in 10% intervals until RH ) 50%, followed by steps of 5% until RH ) 5%. Thereafter, the sorption isotherm was studied using the same steps in RH. Three additional simulations at RH g 40% were started from the final configuration and two intermediate configurations (with stoichiometries of 2.66 and 3.24, i.e., 213 and 259 water molecules in the simulation box) taken from the simulation at RH ) 35%. These simulations were used to specifically probe the stability of the dihydrate form. In addition to the MC moves used for GEMC simulations with fixed unit cell, the simulation with the variable unit cell also employed Parrinello-Rahman moves38,39 to sample individual elements of the H-matrix describing the unit cell shape and volume. For all simulations, the temperature was set to 298 K and the RH is given relative to the saturated vapor pressure for the TIP4P model (psat ) 4.65 kPa),47 which slightly overestimates this quantity compared to experiment. 3. Results and Discussion 3.1. Constant-Stress Monte Carlo Simulations of Anhydrate and Pentahydrate. The unit cell parameters for the anhydrous and pentahydrate forms of compound A obtained from constant-stress Monte Carlo simulations are listed in Table 1. In both cases, the simulations indicate a monoclinic crystal lattice (one angle not equal to 90°), whereas the experiments16 show an orthorombic structure for the anhydrous form and a monoclinic structure for the pentahydrate. The experiments indicate a monoclinic structure for the hemihydrate with a γ angle of 93.48°, i.e., close to that observed in the simulations for the anhydrate. For the pentahydrate, the simulation overestimates the β angle by ≈3°. The mean unsigned errors in the three length parameters are 0.9 and 1.7% for the anhydrate and pentahydrate, respectively. However, the b length is overestimated by about 2% for both forms, whereas the a and c lengths are somewhat underestimated. The unit cell volume is well reproduced for the anhydrous form, but underestimated by 1.5% for the pentahydrate. In order to compare the relative arrangement of compound A molecules in the unit cell between simulation and experiment, the XRD powder patterns were computed by the software package Mercury48 for the ensemble-averaged structures of the anhydrate and pentahydrate. For these calculations the simulated unit cell parameters were adjusted to the corresponding experimental values while the fractional coordinates of the center-ofmass and the orientation for each molecule remained unchanged. The powder patterns (see Figure 3) and overlays of the structures (without scaling of the cell parameters, see Figure 4) indicate a very good match between the simulated and experimental structures; i.e., the placements of the molecules in the unit cells of the anhydrate and pentahydrate are well reproduced by the simulations.


Figure 3. XRD powder patterns measured experimentally15,16 (EXP), calculated for the Rietveld refined structures16 (RV), and calculated for the simulated (SIM) structures of the anhydrous (a) and pentahydrate (b) forms. For clarity, the intensities of powder patterns for the SIM and RV structures are shifted upward by 40 and 20 units, respectively.

3.2. GEMC Simulations with Fixed Unit Cell. A. Moisture Sorption Isotherm. To verify that the simulations with the TraPPE/OPLS force field yield the correct stoichiometry for a given unit cell and to analyze the structures of the hydrates, GEMC simulations with fixed unit cell parameters taken from the experimental structures were carried out. The resulting moisture sorption and desorption isotherms are shown in Figure 5 together with the experimental isotherms.15,16 All four simulated isotherms show typical Langmuir behavior, and the simulations for the hemihydrate, dihydrate, tetrahydrate, and pentahydrate unit cells reach plateau values with stoichiometries of 0.48, 2.3, 4.2, and 5.0 water molecules per compound A molecule at the highest RH. Given the fixed unit cell parameters, the simulations exhibit relatively little hysteresis (only for the dihydrate and tetrahydrate unit cells at RH ) 20%). It is also evident that for all four isotherms the simulations reach a water loading in good agreement with experiment at the RH values where a given hydrate is experimentally observed to be stable; i.e., RH ) 30% for the hemihydrate, 55% for the dihydrate, 70% for the tetrahydrate, and 90% for the pentahydrate. The average numbers of water molecules found in the simulations near these RH values are listed in Table 2 and correspond to stoichiometries of 0.43, 2.1, 4.1, and 5.0 for the four hydrates. It needs to be emphasized that, due to the fixed unit cell parameters used in these simulations, the plateau values are reached at lower RH and persist at higher RH than observed experimentally. B. Structural Analysis. Snapshots of the final configurations are shown in Figure 6. In all cases, the layered structure of the compound A molecules is clearly evident and the relative arrangements of these molecules within the layers are very similar for all hydrates in agreement with the experimental observations. In the tetra- and pentahydrates, the water molecules form two thin layers in the unit cell, and the number of water molecules is similar in each layer of a given hydrate. In the hemihydrate, the neighboring compound A layers are so close

Zhao et al. to each other that the water molecules can only form a channel along the b axis (see Figure 6). The snapshot of the hemihydrate in the ba plane already indicates that water molecules are paired in these channels. The simulated dihydrate structure can be characterized as a hybrid of the hemi- and tetrahydrate structures, i.e., one set of channels and one thin layer leading to unequal occupation between the compound A layers. This structural motif differs from the analysis of the XRD powder pattern that was interpreted as evenly spaced compound A layers which according to a Connolly surface analysis would be separated by incomplete layers of water molecules.15,16 A further set of isotherm simulations started with the dihydrate unit cell parameters and with water arranged between evenly spaced compound A layers, converted throughout the simulations to the hybrid channel/layer arrangement found previously. Thus, this hybrid structure is clearly more stable for the force field used here but not supported by the experimental data. The oxygen-oxygen radial distribution functions (RDFs) for hydration water and the corresponding number integrals (NI) obtained for the four hydrates are depicted in Figure 7. The first peak of the RDFs for the dihydrate, tetrahydrate, and pentahydrate is found at a separation of 2.78 Å, i.e., a distance that falls within 0.03 Å of the first peak position in bulk water at ambient conditions determined experimentally50-52 and obtained from simulations.28,53-55 Due to the thinness of the hydration layers the coordination numbers at rO-O ) 3.5 Å are 2.8, 3.2, and 3.6 for the dihydrate, tetrahydrate, and pentahydrate, respectively, i.e., less than the value for ice-like tetrahedral coordination and the value found for bulk liquid water.28,53-55 For the pentahydrate, the second peak in the oxygen-oxygen RDF is also found close to the bulk-like value of 4.5 Å, whereas the decreasing hydration number leads to a shift of this peak to rO-O ) 4.8 for the dihydrate and tetrahydrate. It should be noted here that the greatly enhanced peak heights compared to bulk liquid water do not reflect a more ordered structure (as can be inferred from the heights of the first minimum that fall above unity for these three hydrates) but are an outcome of the normalization to the overall water densities in the hydrates. The oxygen-oxygen RDF for the hemihydrate differs very significantly from those for the other three hydrates. The first peak is shifted outward to 2.92 Å and the minima are very pronounced, particularly between the second and third peaks and between the third and fourth peaks. The coordination numbers at rO-O ) 3.5, 6.1, and 9.0 Å are only 0.9, 1.8, and 3.5, respectively. The position of the third peak corresponds to the unit cell parameter b. Thus, the water molecules in the hemihydrate strongly interact with only one other water molecule (see below) and their locations are constrained by the compound A molecules. Nevertheless, it should be noted here that the dimer binding energies at rO-O ) 2.78 and 2.92 Å differ by less than 2% for the TIP4P model; i.e., this shift does not lead to a significantly weaker hydrogen bond because the magnitudes of both the favorable Coulomb energy and the unfavorable LJ energy decrease concomitantly. The distributions of the water molecules in the bc plane are shown in Figure 8. The water channels found in the hemihydrate and in one part of the dihydrate stand out very clearly in these projections, and the pairing of water molecules is also evident. Specific sorption sites are less evident for the thin water layers present in the higher hydrates. Nevertheless, these distributions are far from uniform and most peaks reach a height of about 4, i.e., 4 times higher densities than a uniform distribution. The distributions of the azimuthal angle between the dipole vector of water and the three cell axes are presented in Figure



Figure 4. Overlays of the experimental structures15,16 (stick representation) and simulated structures (translucent spheres) of the anhydrous (a) and pentahydrate (b) forms. It should be noted that the unit cell of the anhydrate determined experimentally is orthorombic and contains only four compound A molecules.

Figure 5. Water sorption isotherms obtained from GEMC simulations with fixed unit cells corresponding to the hemihydrate (green), dihydrate (blue), tetrahydrate (magenta), and pentahydrate (orange). The sorption and desorption isotherm data are shown as up and down triangles, and for clarity, only the former are connected by solid lines. Filled and open black squares represent experimental sorption and desorption data.15,16

TABLE 2: Average Numbers of Water Molecules, Nwater, in the Crystalline Simulation Box and of Hydrogen Bonds, B, per Water Molecule with Compound A and with Other Water Molecules, and Isosteric Heat of Adsorptiona hemihydrate dihydrate tetrahydrate pentahydrate % RH Nwater BA Bwater ∆Had [kJ/mol]

40 341 1.191 0.741 -54.15

50 1702 0.471 2.482 -54.22

70 3253 0.371 2.731 -53.33

90 4022 0.271 2.981 -52.72

a Subscripts denote the standard error of the mean in the last digit.

9. The water molecules in the hemihydrate show strong orientational preferences with three maxima for the angles with the a and b cell vectors and two maxima for the remaining angle. In contrast, for the higher hydrates the orientational distributions with the b vector are nearly uniform and those with the a vector show a broad plateau. The distributions with the c vector show two peaks that shift to more perpendicular orientations with increasing hydration number. Overall, it appears that the competition between hydrogen bonding (see below) with the compound A molecules and with neighboring water molecules leads to an orientational disordering of the water molecules for the higher hydrates compared to the hemihydrate.

Further insight into the sorption process can be gleaned from a hydrogen bond (HB) analysis of the hydration water (see Table 2). For the purpose of this analysis, we used a combined distance-angular criterion (rO-O < 3.3 Å and rO-H < 2.5 Å and cos φO-H · · · O < -0.1) that has previously been used for other hydrogen bonding systems.29,56 The total number of HBs per water molecule increases from 1.93 for the hemihydrate to 3.25 for the pentahydrate. Whereas each water molecule in the hemihydrate forms one HB with compound A, only about a quarter of the water molecules form a direct HB to the host for the pentahydrate. This loss in guest-host hydrogen bonding is more than compensated by a quadrupling of the water-water HBs from 0.74 to 2.98. The number of water-water HBs for the hemihydrate is less than onesdespite the pairing of molecules as seen in Figure 8sbecause there are some vacancy defects in the hemihydrate (stoichiometry of 0.85) and because some of the pairs do not satisfy the angular part of the HB criterion. The fact that for the higher hydrates the number of HBs with other water molecules significantly exceeds the number of HBs with compound A may offer an explanation for the less ordered density distribution in the bc plane (see Figure 8). A calculation of the isosteric heat of adsorption for the different hydrate forms indicates that sorption for the hemihydrate is enthalpically slightly more favorable than for the pentahydrate despite the increased number of hydrogen bonds found for the higher hydrate. 3.3. GEMC Simulations with Variable Unit Cell. The simulations with variable unit cell were started from the experimental pentahydrate structure with a run at 100% RH (given the finite size of the simulated system, nucleation does not occur at RH ) 100%). Figures 10 and 11 summarize the outcome of these simulations. It can be seen that, given the freedom to relax the unit cell, the system takes on a little more water and reaches a stoichiometry of 5.2 at RH ) 100% with a concomitant increase in the interlayer spacing, d, by about 1% compared to the experimental pentahydrate structure.15,16 This stoichiometry is also maintained at RH ) 90%. Thereafter, the stoichiometry decreases gradually (and at near constant rate) until a tetrahydrate with a value of 3.9 is reached at RH ) 40%. At his point, the value of d ) 35.1 Å is in excellent agreement with the corresponding experimental value for the tetrahydrate.15 This is followed by an abrupt decrease in the stoichiometry to 1.95 at RH ) 35% and to 0.37 at RH ) 30%, i.e., values corresponding to a dihydrate state and hemihydrate state, respectively. For these two hydrates, the calculated interlayer spacings are 0.7% smaller than the experimental values. Here, it should be noted that these desorption simulations in the RH


Zhao et al.

Figure 6. Snapshots of the final configurations taken from the simulations with fixed unit cell parameters at the RH values given in Table 2. The top row shows the hemihydrate viewed from the a axis (a), the c axis (b), and the b axis (c). The bottom row shows the dihydrate (d), tetrahydrate (e), and pentahydrate (f) viewed from the b axis. For parts b-f, the layered structure is highlighted by aligning the a axis with or close to the vertical direction. The compound A molecules are represented by green wire frames, and water is depicted using red oxygen and white hydrogen atoms. The unit cells are indicated by the white frames.

Figure 7. Oxygen-oxygen radial distribution functions (solid lines) and number integrals (dashed lines) for water in the hemihydrate (black), dihydrate (red), tetrahydrate (green), and pentahydrate (blue) forms at the RH values given in Table 2.

range of 50-30% require extremely long simulations (upward of 106 MC cycles) to reach stable water loadings. This is most likely caused by the need for a concerted decrease in the water loading and the interlayer spacing to extents greatly exceeding the fluctuations encountered frequently throughout the simulations. This points to large free energy barriers49 separating tetrahydrate, dihydrate, and hemihydrate (see below). Once the hemihydrate form is reached, further reduction in RH leads to a gradual decrease in the stoichiometry until a value of 0.03 is reached at RH ) 5%. During this desorption process, the interlayer spacing decreases by another 0.2 Å compared to the hemihydrate found in the simulations at RH ) 30%; i.e., the value of d lies ≈1% below the experimental anhydrate structure. The final structure at RH ) 5% was used as an initial structure for the computation of the sorption isotherm. As can be seen from Figure 10, these simulations track the desorption branch up to RH ) 30%, where a loading of 0.35 is reached. Beyond

this RH, the simulated sorption branch departs from the desorption branch, the loading increases only very slowly, and the hemihydrate form persists up to RH ) 100% with a stoichiometry of 0.48. Thus, the simulations are plagued by a much larger hysteresis loop than observed experimentally.15,16 Clearly, the simulations are not able to surmount the free energy barrier that separates the hemihydrate from the dihydrate. It should be noted here that this transformation would require an increase in the interlayer spacing by about 2.5 Å. In these simulations, the maximum displacement for a cell length move was adjusted to yield a 50% acceptance rate and settled around 0.1 Å, and only single particle exchanges were used. Thus, the simulation cannot “tunnel” through the free energy barrier. Test simulations using (a) larger values for the maximum displacement in the cell parameter, (b) multiparticle swap moves, and (c) concurrent cell parameter and particle swap moves did not lead to improved sampling. Maybe the use of special free energy techniques, such as umbrella sampling57 or metadynamics,58 would help to overcome the sampling problem, but it was not investigated here. To further probe the stability of the dihydrate and tetrahydrate forms that were observed during the desorption simulations, additional simulations were carried out at RH g 40% starting from different initial conditions. First, a set of simulations using the final (dihydrate) structure at RH ) 35% for a sorption run maintains the dihydrate stoichiometry up to RH ) 60%. Even when the RH is further raised to 100%, the stoichiometry increases only to 2.7; i.e., it does not revert to a higher hydrate. Second, simulations at RH ) 40% started with stoichiometries of 2.7 and 3.2 lead to dihydrate and tetrahydrate forms with loadings of 2.3 and 3.8. This indicates that the dihydrate and tetrahydrate are indeed stable and separated by a free energy barrier that has a maximum near a stoichiometry of =3. As also observed in the GEMC simulations with fixed unit cell (see Figure 6), the dihydrate structure found in the simulations with variable unit cell is a hybrid structure containing one layer with water channels (as found for the hemihydrate) and one layer with a thin water film (similar to the higher



Figure 8. Distributions of the fractional center-of-mass location of water molecules in the bc plane for the hemihydrate, dihydrate, tetrahydrate, and pentahydrate (from left to right) at the RH values given in Table 2. The simulated systems consist of five and two unit cells in the b and c directions, respectively. There are two crystallographically independent compound A molecules in the unit cell, there are also two distinct layers and these are represented in the two rows. Please note the difference in the scales for channel and thin-film distributions.

Figure 10. Water sorption isotherms obtained from GEMC simulations with variable unit cell. The desorption and sorption branches are shown as open and filled red circles. The blue and magenta stars indicate additional starting structures used for the sorption simulations depicted by the blue up triangles and the magenta diamonds. The right and left blue triangles indicate the loading of the thin-film layer and the channellike layer found in the hybrid dihydrate structure. Filled and open black squares represent experimental sorption and desorption data.15,16 Figure 9. Distributions of the azimuthal angles between the dipole vector of water and the unit cell vectors a (top), b (middle), and c (bottom) for the hemihydrate (black), dihydrate (red), tetrahydrate (green), and pentahydrate (blue).

hydrates). In the sorption simulations with variable unit cell started from the structure at RH ) 35%, the channel-like layer contains about 1.2 water molecules per unit cell for RH e 60% and 1.6 water molecules at RH ) 100%, whereas the number of water molecules in the thin-film layer increases from 17 to 20; i.e., the latter is close to the stoichiometry found for a single layer of the pentahydrate form. That is, once such a thin-film layer is present, its loading can transform gradually to accommodate loadings equivalent to those found between tetrahydrate and pentahydrate forms. It should be noted that the experimental sorption isotherm also exhibits a smooth transition between these two hydrates and no hysteresis in the RH range from 80 to 90% where the stoichiometry changes from 4.5 to 5.

Furthermore, the gradual decrease in the stoichiometry from the hemihydrate to the anhydrous form (and vice versa) is caused by an increase in the number of lattice defects; i.e., with decreasing RH more and more of the sorption sites seen in Figure 8 are not filled by water molecules. Since this change in the number of occupied sorption sites is not associated with a large change in the unit cell parameters (see Figure 11), there is no hysteresis for RH e 30% and the transformation between anhydrous and hemihydrate forms is gradual. Since all types of experimentally determined hydrate forms were observed in the desorption isotherm, the full range of compound A structures can be compared with experiment. Both simulation and experiment indicate that changes in the water loading result primarily in changes of the unit cell parameter a. A plot of the projection of the a vector onto a vector normal to the bc plane versus the stoichiometry demonstrates that the simulations closely track the experimental behavior for this property (see Figure 11). Linear regressions for the simulation


Figure 11. Projections of the cell vector a onto the normal to the bc plane versus the stoichiometry. Open circles (dashed line) and filled squares (solid line) represent data (and their linear regressions) from the simulated desorption branch and the experimental structures,15,16 respectively.

and experimental data yield y ) 29.52 + 1.41x (with a correlation coefficient of 0.999) and y ) 29.66 + 1.37x (R ) 0.998), respectively; i.e., a unit change in the stoichiometry leads to a change in the interlayer spacing by 1.4 Å. Thus, the simulations reproduce the uniaxial dependence of the unit cell on the stoichiometry (and, hence, the RH) that is most likely responsible for the experimentally observed tablet cracking. 4. Conclusions The structural and moisture sorption properties of the crystalline phases of pharmaceutical compound A have been investigated using Monte Carlo simulations in various ensembles. The development of a formulation for compound A proved problematic because this compound exists in multiple hydrates that lead to tablet cracking at high humidities. Simulations with fixed stoichiometry in the constant-stress ensemble and simulations of the sorption isotherm with fixed unit cell parameters compare well with the experimental data.15,16 That is, the former simulations yield crystalline structures in close agreement with experiment and the latter yield plateau values for the sorption isotherms that possess the correct stoichiometry. Most importantly, simulations of the desorption isotherm in the Gibbs ensemble with variable unit cell are able to predict the existence of all five hydrates (penta-, tetra-, di-, and hemihydrates, and anhydrous form) with structural parameters that satisfactorily match the experiment. However, the simulations suffer from significant hysteresis for RH > 30%. During the desorption branch, the tetrahydrate remains present for a RH as low as 40%, whereas the transformation to the dihydrate occurs already at RH ) 60% for the experimental desorption isotherm.15,16 Furthermore, the desorption simulations yield a dihydrate form only at RH ) 35%, whereas it is found experimentally over the RH range from 55 to 40%. During the sorption simulations, a transformation of the hemihydrate to the dihydrate is not observed even at RH ) 100%; i.e., the simulations are not able to overcome the free energy barrier that separates these structures. As also observed in the experimental work, the transformations from pentahydrate to tetrahydrate and from hemihydrate to anhydrate appear gradually over a range of humidities.

Zhao et al. Structural analysis shows that all hydrate forms have a similar two-dimensional framework consisting of compound A with hydration water in well-defined layers. In the hemihydrate, the water molecules are present as one-dimensional channels oriented along the b axis. However, the mobility of the water molecules in these channels is restricted and the water molecules are found in pairs, i.e., form only one hydrogen bond to another water molecule due to the constraints of the compound A framework. The simulations indicate a hybrid structure for the dihydrate with alternating water channels and thin water films. The water molecules in the latter have considerable lateral and orientational disorder and form at least 5 times as many hydrogen bonds with other water molecules than with compound A. The channel-type water layers are not found for the higher hydrates, and both the tetrahydrate and pentahydrate possess two equivalent layers with thin water films in the unit cell. The structural effect of the changing stoichiometry is dominated by a uniaxial expansion along the a axis of the unit cell with a unity change in stoichiometry resulting in a 1.4 Å change of the interlayer spacing, which provides a molecular level explanation of the tablet cracking. Finally, this study indicates the potential of molecular simulations to be used as a screening tool to probe for the formation of multiple hydrates, in addition to the more common use of polymorph prediction,3-5 before these hydrates may be encountered in later stages of the drug development process. Acknowledgment. We thank Raj Suryanarayanan for helpful discussions. Financial support from the National Science Foundation (CTS-0138393, CBET-0553911, CBET-0756641) is gratefully acknowledged. Part of the computer resources were provided by the Minnesota Supercomputing Institute. Supporting Information Available: Table listing the Lennard-Jones parameters and partial charges used for compound A. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Byrn, S. R. Solid State Chemistry of Drugs; Academic Press: New York, 1982. (2) Bernstein, J. Polymorphism in Molecular Crystals; Oxford University Press: Oxford, 2002. (3) Lommerse, J. P. M. Acta Crystallogr., Sect. B 2000, 56, 697–714. (4) Motherwell, W. D. S. Acta Crystallogr., Sect. B 2002, 58, 647– 661. (5) Day, G. M. Acta Crystallogr., Sect. B 2005, 61, 511–527. (6) Infantes, L.; Fabian, L.; Motherwell, W. D. S. CrystEngComm 2007, 9, 65–71. (7) van Eijck, B. P.; Kroon, J. Acta Crystallogr., Sect. B 2000, 56, 535–542. (8) Chen, J. X.; Wang, J. K.; Zhang, Y.; Wu, H.; Chen, W.; Guo, Z. C. J. Cryst. Growth 2004, 265, 266–273. (9) Wick, C. D.; Siepmann, J. I.; Sheth, A. R.; Grant, D. J. W.; Karaborni, S. Cryst. Growth Des. 2006, 6, 1318–1323. (10) Hulme, A. T.; Price, S. L. J. Chem. Theory Comput. 2007, 3, 1597– 1608. (11) Boek, E. S.; Coveney, P. V.; Skipper, N. T. Langmuir 1995, 11, 4629–4631. (12) Karaborni, S.; Smit, B.; Heidug, W.; Urai, J.; van Oort, E. Science 1996, 271, 1102–1104. (13) Shroll, R. M.; Smith, D. E. J. Chem. Phys. 1999, 111, 9025–9033. (14) Tambach, T. J.; Bolhuis, P. G.; Hensen, E. J. M.; Smit, B. Langmuir 2006, 22, 1223–1234. (15) Kiang, Y.-H.; Xu, W.; Tsou, N.; Ball, R. G.; Stephens, P. W.; Yasuda, N.; DeBusi, L. A. Presented as a poster at PPXRD-3, 2004. (16) Kiang, Y.-H.; Xu, W.; Ball, R. G.; Stephens, P. W.; Yasuda, N. Cryst. Growth Des. 2009, 9, 1833–1843. (17) Cox, J. S. G.; Woodard, G. D.; McCrone, W. C. J. Pharm. Sci. 1971, 60, 1458–1465. (18) Chen, L. R.; Young, V. G., Jr.; Lechuga-Ballesteros, D.; Grant, D. J. W. J. Pharm. Sci. 1999, 88, 1191–1200.

Hydrate Formation in a Pharmaceutical Compound (19) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102, 2569– 2577. (20) Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1999, 103, 4508– 4517. (21) Wick, C. D.; Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 8008–8016. (22) Chen, B.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2001, 105, 3093–3104. (23) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. J. Phys. Chem. B 2004, 108, 17596–17605. (24) Wick, C. D.; Stubbs, J. M.; Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2005, 109, 18974–18982. (25) Lee, J.-S.; Wick, C. D.; Stubbs, J. M.; Siepmann, J. I. Mol. Phys. 2005, 103, 99–104. (26) Jorgensen, W. L.; Briggs, J. M.; Contreras, M. L. J. Phys. Chem. 1990, 94, 1683–1686. (27) Briggs, J. M.; Nguyen, T. B.; Jorgensen, W. L. J. Phys. Chem. 1991, 95, 3315–3322. (28) Jorgensen, J. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926–935. (29) Chen, B.; Siepmann, J. I. J. Phys. Chem. B 2006, 110, 3555–3563. (30) Rafferty, J. L.; Zhang, L.; Siepmann, J. I.; Schure, M. R. Anal. Chem. 2007, 79, 6551–6558. (31) Wick, C. D.; Siepmann, J. I.; Theodorou, D. N. J. Am. Chem. Soc. 2005, 127, 12338–12342. (32) Zhao, X. S.; Chen, B.; Karaborni, S.; Siepmann, J. I. J. Phys. Chem. B 2005, 109, 5368–5374. (33) Wood, W. W.; Parker, F. R. J. Chem. Phys. 1957, 27, 720–733. (34) Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, W. A. Intermolecular Forces; Oxford University Press: Oxford, 1987; p 519. (35) Ewald, P. P. Ann. Phys. 1921, 64, 253–287. (36) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, 1987. (37) McDonald, I. R. Mol. Phys. 1972, 23, 41–58. (38) Parrinello, M.; Rahman, A. Phys. ReV. Lett. 1980, 45, 1196–1199. (39) Yashoneth, S.; Rao, C. N. R. Mol. Phys. 1985, 54, 245–251. (40) Panagiotopoulos, A. Z. Mol. Phys. 1987, 61, 813–826.

J. Phys. Chem. B, Vol. 113, No. 17, 2009 5937 (41) Panagiotopoulos, A. Z. Mol. Phys. 1987, 62, 701–719. (42) Panagiotopoulos, A. Z.; Quirke, N.; Stapleton, M.; Tildesley, D. J. Mol. Phys. 1988, 63, 527–545. (43) Smit, B.; de Smedt, P.; Frenkel, D. Mol. Phys. 1989, 68, 931–950. (44) Siepmann, J. I.; Frenkel, D. Mol. Phys. 1992, 75, 59–70. (45) Mooij, G. C. A. M.; Frenkel, D.; Smit, B. J. Phys.: Condens. Matter 1992, 4, L255-L259. (46) Vlugt, T. J. H.; Martin, M. G.; Smit, B.; Siepmann, J. I.; Krishna, R. Mol. Phys. 1998, 94, 727–733. (47) Chen, B.; Siepmann, J. I.; Klein, M. L. J. Phys. Chem. A 2005, 109, 1137–1145. (48) Macrae, C. F.; Edgington, P. R.; McCabe, P.; Pidcock, E.; Shields, G. P.; Taylor, R.; Towler, M.; van de Streek, J. J. Appl. Crystallogr. 2006, 39, 453–457. (49) Tambach, T. J.; Bolhuis, P. G.; Smit, B. Angew. Chem., Int. Ed. 2004, 43, 2650–2652. (50) Soper, A. K.; Bruni, F.; Ricci, M. A. J. Chem. Phys. 1997, 106, 247–254. (51) Soper, A. K. Chem. Phys. 2000, 258, 121–137. (52) Hura, G.; Russo, D.; Glaeser, R. M.; Head-Gordon, T.; Krack, M.; Parrinello, M. Phys. Chem. Chem. Phys. 2003, 5, 1981–1991. (53) Chen, B.; Xing, J.; Siepmann, J. I. J. Phys. Chem. B 2000, 104, 2391–2401. (54) Kuo, I.-F. W.; Mundy, C. J.; McGrath, M. J.; Siepmann, J. I.; VandeVondele, J.; Sprik, M.; Hutter, J.; Chen, B.; Klein, M. L.; Mohamed, F.; Krack, M.; Parrinello, M. J. Phys. Chem. B 2004, 108, 12990–12998. (55) McGrath, M. J.; Siepmann, J. I.; Kuo, I.-F. W.; Mundy, C. J.; VandeVondele, J.; Hutter, J.; Mohamed, F.; Krack, M. ChemPhysChem 2005, 6, 1894–1901. (56) Stubbs, J. M.; Siepmann, J. I. J. Am. Chem. Soc. 2005, 127, 4722– 4729. (57) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187–199. (58) Martonak, R.; Laio, A.; Parrinello, M. Phys. ReV. Lett. 2003, 90, 075503.

JP808164T

Exploring the Formation of Multiple Layer Hydrates for a Complex

Recommend Documents