Understanding Conformational Polymorphism in Ganciclovir: A Holistic Approach

Abstract

We present a holistic crystallographic study of the antiviral ganciclovir, including insights into its solid-state behavior, which could prove useful during drug development, making the process more sustainable. A newly developed methodology was used incorporating a combination of statistical and thermodynamic approaches, which can be applied to various crystalline materials. We demonstrate how the chemical environment and orientation of a functional group can affect its accessibility for participation in hydrogen bonding. The difference in the nature and strength of intermolecular contacts between the two anhydrous forms, exposed through full interaction maps and Hirshfeld surfaces, leads to the manifestation of conformational polymorphism. Variations in the intramolecular geometry and intermolecular interactions of both forms of ganciclovir were identified as possible predictors for their relative thermodynamic stability. It was shown through energy frameworks how the extensive supramolecular network of contacts in form I causes a higher level of compactness and lower enthalpy relative to form II. The likelihood of the material to exhibit polymorphism was assessed through a hydrogen bond propensity model, which predicted a high probability associated with the formation of other relatively stable forms. However, this model failed to classify the stability of form I appropriately, suggesting that it might not have fully captured the collective impacts which govern polymorphic stability.

1. Introduction

In the last decades, there has been considerable investment into creating sustainable chemical technologies and reactions, specifically within the pharmaceutical industry, in order to maximize the efficiency of the traditionally lengthy and expensive process of drug development [1,2,3]. Computational advancements have been developed in parallel to such efforts, namely to unravel the intrinsic entanglement between the crystalline material and its solid-state properties, both for the possibility of exploiting such knowledge to engineer a pharmaceutical drug with the desired properties, and also to ensure its stability and reliability [3,4]. Given that the orientation and interactions between molecules within a crystal determine a material’s physiochemical characteristics such as solubility, any alterations in this arrangement could have substantial implications, especially in a pharmaceutical context in relation to bioavailability and shelf life. Numerous publications illustrate examples of how different computational techniques are being applied to investigate and predict such features [5,6,7,8,9].

This case study revolves around ganciclovir (Figure 1), which is an acyclic guanine nucleoside analogue, that acts as a DNA synthesis inhibitor to various forms of herpes viruses, thereby reducing the rate of viral growth within the host [10]. In general, ganciclovir is formulated as a sodium salt and is administered through an intravenous route (Cytovene^®) or also as an implant (Vitrasert^®) [11]. It is also used as an ophthalmic gel which contains ganciclovir in its pure form, and a very small percentage of water (Virgan^®) [11,12]. In spite of its importance as a BCS class III antiviral, and most especially as a treatment against cytomegalovirus, ganciclovir’s physicochemical properties have still not been optimized [13,14]. They include poor permeability and limited bioavailability, all of which hinder the drug’s performance. Due to such properties, this antiviral is administered in frequent and high doses, exposing the patient to higher risks of toxicity.

In view of all the above, the need for a complete study arose, entailing a detailed analysis at intramolecular, intermolecular and supramolecular levels using computational methods. This investigation is a comparative study between two anhydrous, enantiotropically related conformational polymorphs I (USP reference standard) and II (Figure 2), in an attempt to identify the factors which influence polymorph stability and assess the possibility of forming other polymorphs [15,16]. Published experimental data characterizing ganciclovir polymorphs are used to complement the observations extracted from the developed methodology [15,17,18], which incorporates a combination of different statistical and thermodynamic approaches, namely the full set of programs and interfaces available through the Cambridge Structural Database (Mercury, Mogul, Isostar) and CrystalExplorer. This approach can be applied not only to the API (Active Pharmaceutical Ingredient) discussed in this article but also to other crystalline materials.

2. Computational Methods

2.1. Non-Bonding Interactions Analysis

The intramolecular geometry analysis was conducted using V1.7.5 Mogul, as described in a paper by Galek et al., together with some adjustments to optimize the results [20]. During this investigation, all organometallic entries were excluded, as their geometric parameters may have introduced a bias in the histograms provided in Mogul when analyzing a purely organic molecule. In cases where few fragments were available, the relevance threshold was reduced from 1.00 to 0.75.

The electrostatic potential was calculated through Mercury by means of the MOPAC (Molecular Orbital Package) interface, using the RM1 semi-empirical Hamiltonian method, and was mapped onto the VdW molecular surface [21,22,23].

IsoStar V2.2.5 was utilized for initial examination of the interactions within ganciclovir, with functional groups acting as the central groups [24]. The analysis of possible intermolecular interactions was taken further by means of the construction of full interaction maps in Mercury, taking into consideration the environmental effects of collective factors and steric exclusion to produce maps that are unique to each form and conformation [25,26].

The detection of any intramolecular and intermolecular hydrogen bonding in the asymmetric unit was investigated using Mercury for both anhydrous forms individually. The definition of a hydrogen bond was modified to have an angle of more than 120°, and the distance range was slightly extended further to be less or equal to the sum of VdW radii + 1.00 Å, to ensure that all potential interactions were included [27].

Hirshfeld surfaces were plotted for each of the anhydrous forms using CrystalExplorer, so as to gain a better understanding of the network of non-bonding contacts, beyond conventional hydrogen bonds [28]. The participants of such contacts were identified through fingerprint scatterplots, which map the distance from a point on the surface to the closest nucleus inside the surface, d_i, against the distance outside the surface, d_e. Further details about the calculation of each descriptive variable used can be found elsewhere [29].

Crystal packing similarity was investigated through Mercury, using both anhydrous forms as reference, and the whole database was explored for any entry having sufficient packing similarity. Default selection options were retained, including a molecular cluster size of 15, 20% distance tolerance and 20° angle tolerance.

2.2. Energy Framework Analysis

Energy frameworks were constructed in CrystalExplorer using the CE-B3LYP model (basis set: 6-31G(d,p)), taking into account only intermolecular interactions within the radius of 3.8 Å from the centroid on the central molecule [30]. All frameworks were scaled equally so as to facilitate comparisons.

2.3. Polymorph Assessment through Hydrogen Bond Propensity Models

Initially, the definition of a hydrogen bond was set as established earlier; however, this resulted in an excessive number of contacts within each form, some of which were chemically incorrect. Therefore, the bond angle was modified to be larger than 133° (the lower limit of the hydrogen bonds detected in both forms). All hydrogen bond acceptors and donors were selected, with the latter list including carbon. Systems with errors or disorder were excluded from the study, together with organometallic compounds and entries with an R-factor > 0.075.

A set of data containing entries whose chemistry was relevant to ganciclovir was generated so as to build the statistical model for prediction. Potential bias or error was decreased by ensuring that each functional group (5 substructures; Please refer to Supplementary Materials) was represented by an adequate number of hits. After analysis, donor or acceptor candidates with a low number of relevant hits were omitted to avoid regression failures. A logistic regression model was fitted on the data, producing its corresponding area under the receiver operating characteristic (ROC) curve value as a measure of the extent of correct predictions. This procedure was performed separately for both polymorphs.

3. Results and Discussion

3.1. Intramolecular Level

In order to initiate the study, Mogul was used to predict the geometric preferences of every bond length, angle, torsion angle and ring within the molecules present in both forms, by accessing a depository of CSD-based libraries [31]. The results for geometric characteristics of form I showed that there were no unusual parameters, as opposed to some properties of form II that were not within the normal CSD distributions. Although these outlying parameters might be a product of poor data quality, which could be inferred from the atomic displacement parameters of form II (Figure 2), some of these unusual geometric characteristics (Please refer to Supplementary Materials) were also observed in other systems involving the ganciclovir molecule, such as the HCl salt and the monohydrate [32].

The electrostatic map of form I (Figure 3) clearly displays the large variety of hydrogen bond participants, with maxima (dark blue) around the hydrogen atoms of the amine group, and minima (dark red) around the oxygen of the carbonyl group and N3 of the imidazole ring, due to the presence of lone pairs of electrons. In addition, the electrostatic map also proves the presence of other potential hydrogen bond contributors, such as the hydroxyl groups, the ether oxygen atom and also carbon atoms. The electrostatic map of form II was very similar to that of form I (Please refer to Supplementary Materials).

3.2. Intermolecular Level

The guanine ring system was one of the substructures identified via IsoStar, illustrated in Figure 4a, with nitrogen and oxygen atoms acting as probes. The distribution of the data points is very close to what was expected, with the majority concentrated around polar contributors. Most of the structures have bifurcated hydrogen bonds, each forming a contact with both O1 and N3. However, the presence of very few turquoise data points indicates that only a minority of hits seems to be at a close distance to the guanine rings (see frequency distribution plot in Supplementary Materials). This observation was unexpected due to the high electronegative character of both nitrogen and oxygen atoms, which in general enables their involvement in relatively strong contacts.

The contour plot in Figure 4b highlights the ability of carbon to act as a hydrogen bond donor, as well as the dominant presence of possible C–H…π interactions or π…π stacking, created due to the delocalized electron density in the aromatic guanine ring system. Other contour plots were constructed using an aliphatic ether as hydrogen bond acceptor and hydroxyl group as the central group. The distributions of the resultant hits of both plots were as expected (Please refer to Supplementary Materials).

Full interaction maps are sensitive to the specific conformation, meaning that each of the forms will have a different map [25]. Comparison of the maps (Figure 5) reveals that the maps associated with form I have a larger area and higher intensity than those pertaining to form II, indicating that the conformation of form I is more accessible for hydrogen bonding. The fully extended conformation of form I enables it to form more hydrogen bonds, which ultimately might be a major contributor to its thermodynamic stability at ambient conditions. Such stability is also evidenced through the fact that form III (hydrate form) converts to form I at temperatures above 180 °C, an exothermal transition which proves the monotropic relationship between the two forms [18]. The lack of accessibility in relation to the ether oxygen in form II, due to the proximity of the hydroxyl group, was earlier highlighted by the orientation of the torsion angle in Mogul results, which clearly demonstrated a degree of folding.

The geometric features of the hydrogen bonds of both forms were compared with those in the literature, and all results were compatible except those involving a C–H donor, which to the best of our knowledge were not included in any published list (Please refer to Supplementary Materials) [17]. The full interaction maps could not predict the involvement of the carbon donors for form I unless the distance levels were significantly increased. Even though the geometric parameters of the contacts with C–H donors are associated with characteristics of weaker interactions, they have a collective effect on the crystal packing and physicochemical properties of form I [33,34]. The addition of methyl carbon and aromatic carbon as separate probes did not alter maps significantly, suggesting that the C–H…π contacts and π…π stacking are less influential in form I, as opposed to in form II (Please refer to Supplementary Materials).

3.3. Hirshfeld Surface Analysis

Hirshfeld surface analysis was performed [9,35,36] in order to gain a better understanding of the differences in hydrogen bond networks and the contribution of weaker contacts within these conformational polymorphs. Since the Hirshfeld surface depends on the spherical atomic electron densities of a particular molecule within its crystal structure, surface differ even between polymorphs [29]. Examination of the Hirshfeld surfaces in Figure 6 shows the distinguishable features of each form, particularly the position and intensity of some of the red areas, which represent close contacts. The Hirshfeld surface of form I is characterized by more red spots relative to that of form II, which as remarked earlier through the full interaction maps, has a folding feature that hinders its accessibility to hydrogen bonds.

The nature of all non-bonding contacts in both polymorphs is presented through the fingerprint plots in Figure 7, by plotting d_i against d_e and hence translating the information provided by the Hirshfeld surface into a 2D format. The plots are overall strongly related, probably due to the similarities between these conformational polymorphs. The greatest proportion of interactions H…H are with a higher percentage contribution (Figure 8) and smaller minimum distance d_i ≈ d_e ≈ 1.1 Å associated with form II. Although such contacts are believed to be repulsive in nature, Matta et al. demonstrated how the net result of their presence has a stabilizing effect of up to a 10 kcal/mol decrease in the total energy of that particular structure [37].

The spikes at the bottom left of both plots are very prominent, exhibiting the dominance of N…H and O…H contacts, due to the highly polar functional groups in ganciclovir. The presence of these contacts is in agreement with the full interaction maps and the hydrogen bonds predicted by Mercury. The greener shades on the spikes of the form I surface plot indicate a higher frequency of contacts having relatively shorter distances. The percentage contribution of O…H interactions is equivalent for both forms, but the global minimum distance d_i + d_e of around 1.7 Å pertains to form I.

The C–H…π contribution is evident in both polymorphs, as shown in the Hirshfeld surface maps; however, it is more significant in form II (Figure 8). Although these interactions are weaker and less directional than the conventional hydrogen bonds, they still contribute to self-assembly and molecular recognition processes, as they reinforce the stability within the supramolecular structure [38]. On the other hand, the absence of large maxima around the area at d_i + d_e ≈ 3.8 Å and of the typical bold red and blue pairs of triangles on the shape index surfaces, indicate the less prominent π…π contacts [39]. Examination of the network of interactions in the two forms, reveals how the side chain prevents close direct stacking of the aromatic systems, hence minimizing the effect of π…π contacts. Any C…C contact seems to be distant and concentrated on the periphery of the guanine ring system, which hints that the aromatic systems are either far apart and/or shifted relative to one another.

3.4. Supramolecular Analysis

Form I is centrosymmetric (P2₁/c) with a unit cell containing two pairs of molecules that are identical within the pair and inversely-related between pairs. Form II is non-centrosymmetric (P2₁), thus lacking an inversion center.

The extensive hydrogen bonding network in form I is highly visible along its packing pattern, which seems to follow zigzag lines when viewed along c (Figure 9a). This type of pattern allows the layers of molecules to be very close to one another, hence enabling the optimal use of space available [40]. While the molecular volume of form I is greater than that of form II, packing in the former is more compact as evidenced by the volume of its Hirshfeld surface, packing coefficient and density value (

Table 1. Numerical Hirshfeld surface packing data for forms I and II.

	Form I	Form II
Volume (Å3)	250.18	268.07
Packing coefficient	0.770978	0.725452
Density (g/cm3)	1.653	1.549

). Figure 9b shows where the purine backbones seem to be parallel to one another, while the side chains are “out of plane”. Layering of the guanine rings creates an off-center parallel stacking, but distance and angular parameters of the guanine ring centroids are at the upper limits for effective π…π interactions, suggested by the literature [41,42]. This property was also highlighted by a low percentage contribution of C…C contacts in the fingerprint plots.

The packing arrangements in both forms were compared to other crystal structures with similar molecular arrangements. However, even though there were five more entries in the CSD which involved ganciclovir as a component, only the two hydrates (mono- and tri-), resulted in having one out of a cluster of 15 molecules, with a degree of similarity when compared to both anhydrates [32] demonstrating the uniqueness of the packing patterns of form I and II. This test also verified the significant differences between the way in which both anhydrous forms were arranged three-dimensionally resulting in an RMS of 1.017 (see Figure 8).

3.5. Energy Frameworks

The thermodynamic stability profile of each of the polymorphs is a direct implication of the nature and strength of their non-bonding interactions. Energy frameworks were calculated in an attempt to understand the topological dissimilarities of the energy components of the two polymorphs, and subsequently potentially link these characteristics to their respective packing and thermal behavior [35]. The total energy is divided into Coulomb forces/electrostatic potential forces and dispersion energy.

The negative energy components in form I are distributed along both parallel and non-parallel molecules (Figure 10, top row), hence stabilizing the energy architecture throughout the supramolecular structure in multiple directions. The dispersion forces can mainly be traced along the off-centered parallel stacking (along the axis a), as well as between adjacent molecules whose aromatic systems lie along the same plane. The strongest energy contribution (−164 kJ/mol) is attributed to N1–H…O2 and N2–H…O4 interactions, located between parallel molecules. The strength of these stabilizing interactions is due to both electrostatic and dispersion forces, the former being the dominant one. The O4–H…O3 contact contributes towards a lowering of energy (−56.4 kJ/mol), while the combination of O3–H…N3 and C6–H…O1 interactions induces a further stabilizing effect of −63.7 kJ/mol.

In form II, N1–H…O1 and N2–H…N3 are the strongest interactions both of which contribute towards the dominant electrostatic energy (−73.6 kJ/mol) along the planes of aromatic rings. A secondary energy contribution (−45.2 kJ/mol) in the same direction is attributable to the N1–H…O4 contact. The out-of-plane interactions are mainly a result of dispersion contacts and the O4–H…O1 electrostatic interaction. However, such a contribution is considerably low (−18.3 kJ/mol) relative to the other energy components along the aromatic systems plane.

The less exhaustive network of intermolecular bonding in form II, in conjunction with the presence of an intramolecular bond, seems to enhance stability within the molecule. In contrast, form I entails a very complex extensive framework of intermolecular contacts, with 16 hydrogen bonds per molecule (and eight different neighboring molecules). Such a multidirectional framework creates a greater stabilizing effect in form I, as illustrated by the relatively thicker cylindrical radius of the total energy tubes (Figure 10). Such an influence, coupled with a more compact packing arrangement, might contribute towards a lower enthalpy in form I relative to form II, which according to L. Yu must be accompanied by a lower entropy in order to conserve the enantiotropic relationship between the two polymorphs [18,43]. With an increase in temperature, in the range of 222 °C to 228 °C, enough energy is absorbed to break the intermolecular framework in form I, which then transitions to form II (as confirmed by variable temperature X-ray diffraction and differential scanning calorimetry) [15,17,18].

The relatively loose packing and the intermolecular framework in form II seems to be constructed in such a way as to cater for higher energy environments. The energy framework of this polymorph reveals how the stabilizing effect is expanded along planes, whereas much fewer interactions are observed between layers. This lower extent of intermolecular forces between parallel planes in form II might create an environment that can accommodate higher entropy within the crystal structure, thereby allowing it to be thermodynamically stable at higher temperatures.

3.6. Polymorph Assessment

The objective of this analysis revolved around stability assessment of the hydrogen bond network present in both forms, given the complexity of the chemical environment of ganciclovir. The reliability and predictive ability of the fitted logistic regression model were mirrored by the value of the area under the ROC curve, which was equal to 0.878, and by the reduction from the null to the residual deviance [44,45]. The same procedure was performed for each form individually, and since they are polymorphs of the anhydrous form of ganciclovir, very similar coefficients were obtained (Please refer to Supplementary Materials). Minor differences can be attributed to the non-deterministic nature of the process involving the fitting of the model. This outcome confirms the robustness of this method, which is capable of assessing the stability of various polymorphs having different hydrogen bonds simultaneously.

The final model was used to calculate the propensities of all possible intermolecular hydrogen bonds, with those in form II having the highest probabilities (more conventional contacts) (Please refer to Supplementary Materials). The overall low likelihood of the intermolecular bonds in form I was illustrated in the putative structure landscape, which categorized this polymorph as having the least stable hypothetical forms (Please refer to Supplementary Materials). The utilization of every functional group in form I, for intermolecular bonding, might have been prioritized over the formation of the fewer and more probable contacts.

It is understandable that due to low frequencies, the contacts in form I were ranked as having lower probability, and hence low stability was predicted. However, this outcome was not in agreement with the fact that form I is the thermodynamically stable polymorph, under ambient conditions. At this point, it is essential to recall the statistical mechanism behind the construction of the Hydrogen Bond Propensity (HBP) model, which is based on the occurrence of hydrogen bonds in similar chemical environments. Therefore, in some examples, the resultant model might not be sufficient to explain the complexity of hydrogen bonding and to capture the collective effects of multiple factors that determine the polymorph stability [6]. In his research paper, Abramov commented how in general, a HBP model cannot account for enantiotropic relationships between polymorphs, such as the one between form I and II [6]. Moreover, one has to take into account that the directional features and geometric parameters of the contacts are not being considered in the model, and these characteristics have a significant effect on stability.

In an attempt to overcome these limitations, the HBP model was constructed using a much larger training set, but similar results were obtained. Despite such limitations, there was still valuable information that could be extracted from these results. The putative structure landscape in Figure 11 portrays the presence of data points located very close to form II, representing reasonable hypothetical structures having very strong hydrogen bond interactions. The viability of the formation of such forms is highly encouraging in view of further research dedicated to the exploration of other possible ganciclovir polymorphs.

4. Conclusions

This investigation presents an extensive study of ganciclovir from a crystallographic point of view, ranging from a structural to energy framework analysis. The methodology developed has proven how the use of multiple computational tools which target the examination of different characteristics, can be successful in providing a reasonable insight into the mechanism of polymorphism. Furthermore, it can be implemented to extract information about the influence of specific bonds and substituents on the formation of a particular form, regardless of the complexity or nature of the molecule.

Such a dominant impact was observed through the folding of the chain moiety in form II, which intrinsically affects its intramolecular and intermolecular characteristics. On the other hand, full interaction maps and Hirshfeld surfaces demonstrated how the extended chain moiety in form I enables the molecule to have an extensive network of non-covalent contacts, with a significant degree of directionality. These properties were identified as significant factors responsible for the known thermodynamic behavior of both forms, ultimately resulting in the conformational polymorphism present between them. The collective analysis of results from all molecular levels and energy frameworks provides a plausible explanation for the stability of both polymorphs at different temperature ranges.

The nature of this methodology makes it applicable to a wide range of crystal forms, even beyond active pharmaceutical ingredients. Such knowledge would facilitate the selection of co-formers or environmental conditions to selectively target the formation of a particular co-crystal. It would also aid navigation through the possibility of forming other relatively stable polymorphs, as suggested by the putative landscape in the case of ganciclovir. Therefore, the utilization of such data could lead to a more sustainable process of drug development, as well as possibly improving the pharmacokinetic properties of this active pharmaceutical ingredient.