The global energy minimum of DDS, obtained using Gaussian09 (Frisch et al., 2009), was used in the CSP searches. 350,000 Z′ = 1 monohydrate structures were generated using CrystalPredictor2.0 (Karamertzanis and Pantelides, 2005, 2007; Habgood et al., 2015) in 48 common space groups for organic molecules (Supplementary Material). The molecules were held rigid and the lattice energy was evaluated by an exp-6 potential with atomic charges derived using the CHELPG scheme (Breneman and Wiberg, 1990) and minimized. The 10,000 lowest energy crystal structures were used as starting points for optimizing the intermolecular lattice energy (U_inter), with an improved model for the intermolecular forces. This was calculated using the FIT exp-6 potential parameters (Coombes et al., 1996), the sulfur potential derived by Scheraga (Day et al., 2009) and the distributed multipoles (Stone, 2005) derived from the PBE0/6-31G(d,p) charge density using GDMA2 (Stone, 2010).

The optimal proton positions of the amino group and orientation of the phenyl groups, in all crystal structures within 15 kJ mol⁻¹ of the global minimum (116 structures), were determined using the CrystalOptimizer database method (Kazantsev et al., 2011). This was done by minimizing the lattice energy (E_latt), calculated as the sum of the intermolecular contributions (U_inter) and the conformational energy penalty paid for distortion of the molecular geometry to improve the hydrogen bonding geometries. Conformational energy penalties (ΔE_intra, with respect to the pyramidal global conformational energy minimum) and isolated molecule charge densities were computed at the PBE0/6-31G(d,p) level, for each conformation considered in the minimization of E_latt. All isolated-molecule wave function calculations were performed using Gaussian09 (Frisch et al., 2009) and intermolecular lattice energies using DMACRYS (Price et al., 2010).

The 100 most stable structures (within 30 kJ mol⁻¹ of the global minimum) were used as starting points for periodic electronic structure calculations. The DFT-D calculations were carried out with the CASTEP plane wave code (Clark et al., 2005) using the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) exchange-correlation density functional (Perdew et al., 1996) and ultrasoft pseudopotentials (Vanderbilt, 1990), with the addition of a semi-empirical dispersion correction developed by Tkatchenko and Scheffler (TS) (Tkatchenko Scheffler and Scheffler, 2009). Brillouin zone integrations were performed on a symmetrized Monkhorst–Pack k-point grid with the number of k-points chosen to provide a maximum spacing of 0.07 Å⁻¹ and a basis set cut-off of 780 eV. The self-consistent field convergence on total energy was set to 1 × 10⁻⁵ eV. Energy minimizations were performed using the Broyden–Fletcher–Goldfarb–Shanno optimization scheme within the space group constraints. The optimizations were considered complete when energies were converged to better than 2 × 10⁻⁵ eV per atom, atomic displacements converged to 1 × 10⁻³ Å, maximum forces to 5 × 10⁻² eV Å⁻¹, and maximum stresses were converged to 1 × 10⁻¹ GPa. Isolated molecule minimizations to compute the isolated DDS and water energy (U_gas) were performed by placing a single molecule in a fixed cubic 35 × 35 × 35 ų unit cell and optimized with the same settings as used for the crystal calculations.

The experimental hydrate, lower hydrates and forms II and III of DDS, as well as other selected hydrates (Supplementary Material), were minimized with CASTEP and the same settings were used as described for generating the monohydrate crystal energy landscape.

The pair-wise energy contributions to 0.33-Hy and other well-characterized hydrate structures, have been calculated using CrystalExplorer V17 (Turner et al., 2014, 2015; Mackenzie et al., 2017). The optimized atomic positions (PBE-TS) have been used in all subsequent intermolecular interaction energy calculations. The model energies have been calculated between all unique nearest neighbor molecular pairs. The used model (termed CE-B3LYP) uses B3LYP/6-31G(d,p) molecular wave functions calculated by applying the molecular geometries extracted from the crystal structures. This approach uses electron densities of unperturbed monomers to obtain four separate energy components: electrostatic (E_E), polarization (E_P), dispersion (E_D), and exchange-repulsion (E_R). Each energy term was scaled independently to fit a large training set of B3LYP-D2/6-31G(d,p) counterpoise-corrected energies from both organic and inorganic crystals. The CE-B3LYP energies reproduced the training set energies with a mean absolute deviation of ~1 kJ mol⁻¹ (Turner et al., 2014).

Conformational energy scans were performed at the B3LYP/6-31G(d,p) level of theory using Gaussian09 (Frisch et al., 2009), allowing the two torsion angles defining the position of the phenyl rings, C–C–S–C, to rotate by 360° in 20° steps.

Dapsone form III (purity 97%) was purchased from Aldrich. The obtained sample was recrystallized from a hot-saturated methanol solution. The solid product was isolated by filtration and consisted of form III. The organic solvents used were all of analytical grade and purchased from Aldrich or Fluka.

DDS 0.33-Hy was prepared as follows: (i) a slurry of DDS form III in water was stirred in the temperature range from 10 to 30°C for 1 week. The suspension was filtered and the solid was stored at ambient conditions. (ii) A hot saturated solution of form III in water (close to the boiling point) was cooled to room temperature (RT). Within 2 days large elongated 0.33-Hy crystals were obtained.

For hot-stage thermomicroscopic (HSM) investigations a Reichert Thermovar polarization microscope, equipped with a Kofler hot-stage (Reichert, A), was used. Photographs were taken with an Olympus DP71 digital camera (Olympus, A).

Differential Scanning Calorimetry (DSC) thermograms were recorded on a Diamond DSC (Perkin-Elmer Norwalk, Ct., USA) controlled by the Pyris 7.0 software. Using a UM3 ultramicrobalance (Mettler, Greifensee, CH), samples of ~5–7 mg were weighed into perforated or sealed aluminum pans. The samples were heated using rates in between 1 and 20°C min⁻¹ and cooled using a rate of 5 or 10°C min⁻¹ with dry nitrogen as the purge gas (purge: 20 mL min⁻¹). The instrument was calibrated for temperature with pure benzophenone (mp 48.0°C) and caffeine (236.2°C), and the energy calibration was performed with indium (mp 156.6°C, heat of fusion 28.45 J g⁻¹). The errors on the stated temperatures (extrapolated onset temperatures) and enthalpy values were calculated at the 95% confidence interval (CI) and are based on at least five measurements.

Thermogravimetric Analysis (TGA) was carried out with a TGA7 system (Perkin-Elmer, Norwalk, CT, USA) using the Pyris 2.0 Software. Approximately 7–10 mg of sample was weighed into a platinum pan. Two-point calibration of the temperature was performed with ferromagnetic materials (Alumel and Ni, Curie-point standards, Perkin-Elmer). Heating rates of 5 and 10°C min⁻¹ were applied and dry nitrogen was used as a purge gas (sample purge: 20 mL min⁻¹, balance purge: 40 mL min⁻¹).

PXRD patterns were obtained using an X'Pert PRO diffractometer (PANalytical, Almelo, NL) equipped with a θ/θ coupled goniometer in transmission geometry, programmable XYZ stage with well plate holder, Cu-Kα_{1, 2} radiation source with a focusing mirror, a 0.5° divergence slit, a 0.02° Soller slit collimator on the incident beam side, a 2 mm antiscattering slit, a 0.02° Soller slit collimator on the diffracted beam side and a solid state PIXcel detector. The patterns were recorded at a tube voltage of 40 kV and tube current of 40 mA, applying a step size of 2θ = 0.013° with 200 s per step in the 2θ range between 2° and 40°. For non-ambient RH measurements, a VGI stage (VGI 2000M, Middlesex, UK) was used.

The diffraction patterns were indexed using the first 20 peaks with DICVOL04 and the space group was determined based on a statistical assessment of systematic absences (Markvardsen et al., 2001) as implemented in the DASH structure solution package (David et al., 2006). Pawley fits (Pawley, 1981) and Rietveld refinements (Rietveld, 1969) were performed with Topas Academic V5 (Coelho, 2012). The background was modeled with Chebyshev polynomials and the modified Thompson-Cox-Hastings pseudo-Voigt function was used for peak shape fitting. For the Rietveld refinements the DDS and water molecules were treated as rigid body molecules using the PBE-TS optimized conformations of the 0.33-Hy structure.

Moisture sorption and desorption studies were performed with the automatic multisample gravimetric moisture sorption analyser SPS23-10μ (ProUmid, Ulm, D). Approximately 500–750 mg of sample was used for each analysis. The measurement cycles were started at 60% with an initial stepwise desorption (decreasing humidity) to 0%, followed by a sorption cycle (increasing humidity) up to 90% relative humidity (RH), a desorption cycle to 0% RH and a final sorption cycle to 90% RH. The RH changes were set to 2% and the equilibrium condition for each step was set to a mass constancy of ± 0.001% over 60 min and a maximum time limit of 48 h per step.

DDS form III was stirred (500 r.p.m.) in 1.5–2.5 mL of each methanol and water mixture [each containing a different mole fraction of water corresponding to a defined water activity Zhu et al., 1996, Supplementary Material] at 25.0 ± 0.1°C for 21 days. Samples were withdrawn, filtered and the resulting phase was determined using PXRD.

The fact that the asymmetric unit of 0.33-Hy (Kuz'mina et al., 1981; Bel'skii et al., 1983; Yathirajan et al., 2014) consists of four crystallographically distinct molecules (three DDS and one water) makes CSP studies for the experimental stoichiometry too time-consuming. Therefore, we decided to generate the monohydrate crystal energy landscape (one DDS and one water molecule) with the aim to estimate whether water molecules can compete against the DDS-DDS intermolecular interactions and form strong DDS-water contacts. In Figure 2 the computed monohydrate structures are plotted according to lattice energy, which equals the energy that would be required to separate the molecules to infinity, against packing index.

To estimate whether any of the hypothetical monohydrate structure is competitive in energy with form III, we compared the lattice energy of the hydrate (E_latt-Hy) to the lattice energies of the anhydrate (E_latt-III) and ice (E_latt-ICE). If E_latt-Hy < E_latt-III + E_latt-ICE (we assume that hydrate formation is thermodynamically driven), then the hydrate is more stable than the anhydrate. Using the lattice energies of 0.33-Hy, form III (Table 1) and a value of −59 kJ mol⁻¹ (Whalley, 1957, 1976) for ice, as the used functional is known to overbind the ice crystal structures (Thierfelder et al., 2006; Beran and Nanda, 2010), then only one structure, 01_1963, was calculated to be more stable than III. The most stable hypothetical monohydrate was estimated to be 12.38 kJ mol⁻¹ more stable than form III, which is a respectable potential energy differences (Δ_trsU) for a monohydrate with respect to an anhydrate. Thus, the CSP study clearly indicates hydrate formation.

Furthermore, Δ_trsU for the 0.33-Hy to form III was calculated to be 15.37 kJ mol⁻¹, indicating that the experimental hydrate is 3 kJ mol⁻¹ more stable than the computed lowest energy monohydrate structure and rationalizing why the stable 0.33-hydrate and not a monohydrate is formed experimentally.

To confirm that 0.33-Hy is indeed the stable DDS hydrate form and that 01_1963 is not a yet undiscovered monohydrate we subjected DDS to an experimental hydrate screening program. Evaporative crystallization experiments of DDS from a saturated (20°C) aqueous solution, as well as cooling crystallization experiments from hot (boiling) saturated solutions in water at 5°, 25°, 50°, and 75°C resulted in 0.33-Hy crystals in the form of elongated plates. In contrast, evaporation experiments of a hot-saturated solution of DDS in water resulted in a mixture of 0.33-Hy and form III. Slurry experiments in water, isothermal or cycling between 5° and 50°C, always yielded the 0.33-Hy.

Another successful way to produce hydrates are moisture sorption experiments. Therefore, form III and V (see section Moisture Dependent Stability of the Hydrate) were subjected to automated and manual water vapor sorption experiments. Neither form III, nor form V showed a transformation in the RH range up to 90%. Furthermore, no transformation was observed in long-time storage experiments of the two anhydrous DDS forms over saturated KOAc (24% RH), K₂CO₃ (43% RH), NaCl (75% RH), KNO₃ solutions (92% RH) or water (100% RH) within 3 months (end of experiments) at RT and 8°C. Similarly, also 0.33-Hy did not dehydrate or transform to another hydrate if stored under the same conditions over the same time period.

To understand the nature and stability of a hydrate it is important to consider the location and interaction energies of the water molecules in the framework of the host structure. If water molecules are located in open structural voids (tunnels or connected pockets) the term channel hydrate (Brittain et al., 2009) is commonly used. In such hydrates the water molecules may be mobile and may readily escape through these tunnels on modest increase in temperature or decrease in relative humidity (RH). In contrast, if the water molecules are located at isolated-sites (Brittain et al., 2009), it is assumed that water egress is not as facile and requires a considerable rearrangement of the hydrate packing to allow the release of water molecules. This rearrangement results mostly in the formation of a different packing arrangement or in a partial or total collapse of the structure yielding a disordered or amorphous state upon dehydration. The non-stoichiometric behavior of hydrates, where the water content in the structure depends on the water vapor pressure of the surrounding medium (atmosphere), is normally observed in channel hydrates and not in isolated-site hydrates. However, as demonstrated above, DDS 0.33-Hy shows clearly the typical features of an isolated-site hydrate (Figure 3A) but on the other hand shows a non-stoichiometric (de)hydration behavior (see Figure 7) which is a contradiction and questions the common relation between structural features and the stability of hydrates.

To further clarify why the hydrate water can escape easily from the “isolated sites” in the 0.33-Hy, without disrupting the structure, we estimated the pair-wise interaction energies for DDS and water molecules in the 0.33-Hy structure (Table 2). The use of the CE-B3LYP energies and not a specific water potential, not considering specific effects (nanoscale dielectric responses of water) and three-body energy terms were justified by the fact that the contribution of the water to 0.33-Hy was found to be in reasonable agreement with the experiment (Table 3). The water interactions contribute −13.85 kJ mol⁻¹ to the 0.33-Hy lattice. The CE-B3LYP energies of 0.33-Hy (−142.87 kJ mol⁻¹) and Hy_dehy (−132.00 kJ mol⁻¹, optimized RT structures) differ by 10.87 kJ mol⁻¹, ignoring the conformational changes which are expected to account for < 0.2 kJ mol⁻¹ in the case of DDS. Thus, the sum of the interaction energies (E_cluster) roughly corresponds to the lattice energy. Furthermore, we calculated the intermolecular energies for water-host, water-water and host-host molecules for a series of well characterized organic hydrate systems (pharmaceuticals and model compounds) and contrasted the values to the DDS 0.33-Hy (Figure 11, Supplementary Material). The chosen test set consists of stoichiometric and non-stoichiometric hydrates, as well as of channel and isolated-site hydrates. This analysis should indicate, whether it is possible to assess hydrate stability and/or the dehydration mechanism from general features of the hydrate structure and not the location of the water molecules in the structure alone.

In Figure 11 the calculated interaction energies are grouped into contributions arising from host-host interactions and the incteractions including water molecules (host-water and water-water). Furthermore, the ratio between the host-host and water interactions has been calculated. The hydrates are ranked according to the compound:water ratio. For three of the chosen hydrates (dapsone, indinavir, and brucine) it is possible to remove the water molecules under maintaining the crystal structure, which is characteristic for a non-stoichiometric dehydration behavior. A grouping into channel and isolated-site hydrates is not always straight forward, in particular if more than one water molecule is present (di-, tri-hydrate, etc.,). However, based on the energetic contributions, host-water vs. water-water, such a classification is facilitated as it can be expected that in isolated-site hydrates the host-water and in channel hydrates the water-water interactions predominate, respectively. Furthermore, the sum of host-host interactions are stronger in isolated-site hydrates than in channel hydrates.

A requirement for maintaining the crystal lattice upon dehydration is that the hydrate structure exhibits strong and/or a predominance of host-host interactions. Indeed, the three non-stoichiometric hydrates of the test set show the highest percentage of host-host interactions, i.e., ~90% of the interaction energies for the DDS and indinavir hydrates. This value is lower for brucine (61.5%) but compared to dihydrates showing a stoichiometric behavior, brucine exhibits the most/strongest host-host interactions. The fact that the DDS and indinavir hydrates are isolated-site hydrates and brucine is a channel hydrate highlights that it is not possible to deduce whether a stoichiometic or non-stoichiometric dehydration mechanism occurs from the location of the water molecules in the hydrate structure alone. Though, it is possible to rationalize a non-stoichiometric dehydration behavior from the energy contributions of the intermolecular interactions considering the compound:water ratios. On the other hand, the analysis shows that the moisture- or temperature dependent stability of hydrates cannot be derived from the interaction energy calculations. For example, the 5-flucytosine monohydrate (I) already dehydrates at RH values < 40%, whereas 4-aminoquinaldine monohydrate (Hy1A) dehydrates only at RH values below 10% (RT), but exhibits less energetic contributions from the host-host interactions than the 5-flucytosine monohydrate (I).

In the case of 0.33-Hy the latter analysis (Figure 10) strongly indicates that water is only weakly bound and rationalizes the facile moisture- and temperature dependent dehydration behavior.

Knowledge of how water vapor is sorbed by a hygroscopic material and how moisture affects the physical and chemical stability of a (pharmaceutical) product is a crucial question in developing drug products or preparations produced from other fine chemicals. Failures and time delays in product developments can be minimized or avoided with knowledge compiled in thorough solid state investigations. Hydrates require a thorough evaluation of their composition and stability under production relevant conditions and additionally the transformation pathways between different solid state forms of a compound, as well as their stability ranges, should be elucidated. This is mandatory to select the ideal solid state form that guarantees an optimal product performance and stability. In general, non-stoichiometric hydrates are undesired solid forms because any change in water vapor pressure of the surrounding medium causes a change in the water content of the substance, which can be critical for weighing and dosing operations and may thus lead to errors in any analyses, which require exact sample amounts. Such variations in the water content are often difficult to avoid as it requires special efforts to precisely control temperature and humidity conditions during processing and storage. Furthermore, the water molecules which have been released from such a hydrate may interact with other excipients in a drug formulation. Gravimetric moisture sorption/desorption studies (Figure 7), combined with environmental PXRD experiments (Supplementary Material) are the preferred analytical techniques for unraveling this non-stoichiometric behavior of a hydrate.

DDS 0.33-Hy is a prime example for an isolated site hydrate with non-stoichiometric dehydration behavior. The latter behavior may be expected for hydrates where the water is located in open voids such as channels or layers. Thus, this study highlights that the popular structural classification of hydrates into isolated-site hydrates (water molecules are isolated from direct contact), channel hydrates (chains of water molecules) and ion-associated hydrates (metal ions are coordinated with water) cannot be directly related to the dehydration behavior or dehydration mechanism of a hydrate. However, by complementing the structural features with intermolecular energy calculations the observed dehydration behavior can be rationalized. As shown for DDS 0.33-Hy the water molecules are only weakly bound (Figure 11), allowing a facile water egress/ingress with changing environmental conditions. Furthermore, the energy difference between the isomorphic dehydrate structure and anhydrate polymorphs is small. The lattice energy difference between Hy_dehy and form II was calculated as 1.6 kJ mol⁻¹ (PBE-TS, at −273°C) and the transition enthalpy between Hy_dehy and form II determined to be 1.08 ± 0.05 kJ mol⁻¹ (experimentally measured at ~100°C). Thus, the calculations rationalize and indicate the non-stoichiometric dehydration mechanism.

Modeling and predicting hydrate structures of pharmaceuticals still represent a big challenge in computational chemistry. Numerous potentials have been developed for modeling water (Guillot, 2002), however, there exists no method that can sufficiently model all its abnormalities. In organic hydrates the water molecules may be described as confined at nanoscales, implying frustration in their hydrogen-bonding coordination. Consequently, accurately modeling the water molecules in a hydrate lattice requires modeling efforts which go well beyond the methods applied for modeling the organic solid state (Reilly et al., 2016). In lead optimization is has been demonstrated that incorporating the three-body energy terms, and modeling frustration and frustration-related dielectric responses, significantly improves the results (Fernández, 2016, 2017; Fernandez and Scott, 2017). Considering and modeling the latter can be expected to significantly increase the accuracy of lattice and intermolecular energy calculations of water containing species, albeit at the expense of computational cost.

A major difficulty with using CSP in hydrate solid form screening is the computational expense in time and resources to generate the crystal energy landscape for all possible hydrate stoichiometries (mono-, di-, etc.,). However, solid form modeling at the electronic and atomistic level can provide vital support for unraveling the solid state for a compound which may not be achieved with experiments alone. A CSP study answers the question what types of crystal packings are favorable for a specific molecule, unraveling the compromises between close packing efficiency, conformational preferences and the different types of intermolecular interactions that can lead to feasible structures for a molecule (polymorph) or multi-component system (salt, solvate, hydrate, co-crystal). It should be stressed that CSP aids the interpretation of the experimental data (Price et al., 2014) and can guide experimentalists to find new solid forms (Arlin et al., 2011; Braun et al., 2014b, 2016; Neumann et al., 2015; Srirambhatla et al., 2016).

To significantly reduce the computational cost, and to make the calculations feasible, we did not attempt to computationally screen for different hydrate stoichiometries (Braun et al., 2011b) for the chosen model compound DDS, but used the crystal energy landscape of the 1:1 stoichiometry (monohydrate) as a guidance for hydrate formation. The monohydrate crystal energy landscape (Figure 2) shows only one hydrate structure that is more stable than the non-solvated form III and thus indicates hydrate formation.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.