Theoretical search of crystal polymorphs of temozolomide.

Possible polymorphic forms of the chemotherapy drug, temozolomide were predicted from the ab initio and DFT methods. The lattice minimization via distributed multipole analysis was carried out for the hypothetical generated structures. A crystal with unit cell parameters close to the real one and of same space group was retrieved, with partly similar packing and interactions. The analysis of inter molecular interaction (through Hirshfeld surface) and electrostatic potential reveals the complementary sites in the molecule. The 26 predicted structures were analyzed with respect to two computed lattice energies and hydrogen-bond propensity. The lattice energy of the real crystal [EXP] packing ranked number 6 compared on the basis of DMACRYS software and number 3 on the basis of the total lattice energy issued from the Crystalexplorer17 software at the B3LYP/6-31G∗∗ level of theory. The molecule has two strong hydrogen bond donors and five strong acceptors. The predicted packings are stabilized by one or two strong N-H…O/N-H…N as well as weak C-H…O/C-H…N and H…π hydrogen bonds. While the real structure with Z' = 1, EXP, forms only one strong H-bond (N-H…O=C), several of the predicted packings form two strong H-bonds. Two predicted crystal packings have unit cell parameters close to the real structure, one of them shares several common intermolecular interactions.

Introduction

The phenomena of Polymorphism has acquired great academic and industrial interest [1] and the physical and chemical properties of polymorphic structures make a significant impact in both experimental and theoretical investigations. The polymorphic behavior in pharmaceutical compounds leads to different bio-availabity. The theoretical search for unknown polymorphic structures has gained its importance as their experimental confirmation for their existence is the difficult procedure [2]. In this view, the current study aims to explore the different crystal phases of temozolomide drug which are thermodynamically stable, along with the experimentally known polymorphs. For the treatment of malignant glioma, the United States Food and Drug and Administration (US FDA) [3] approved the drug Temozolomide, an oral alkylating agent. The structure of temozolomide molecule [Figure 1] has tetrazine and imidazole rings, labelled as 1 and 2 with ketone and amide functional groups, which have the potential to form hydrogen bonds. The N–H…N and N–H…O hydrogen bonds are the lead interactions for the polymorphic formation. Ten polymorphs of TMZ were reported in the patent literature [4, 5]. Lowe et al. (2008) reported a first crystal structure of temozolomide (called EXPb) with Z’ = 2. Babu et al [5] reported the same crystal form, one polymorph with Z’ = 1 (EXP) and one polymorph with Z’ = 2 (EXPb). And the protonated form of TMZ was also reported [6] in which the crystal structure was solved with one neutral temozolomide, one temozolomide-H+-Cl-, one H3O+.Cl− and three water molecules. They described also three cocrystals of temozolomide. In addition, temozolomide polymorphs, characterized only by powder diffraction, were disclosed in two patents [7] and several of these polymorphs lack crystal structural information. To gain structural insights, various polymorphic forms of temozolomide were predicted theoretically via ab initio calculations [8].

Figure 1

Chemical structure of Temozolomide.

Materials and methods

Crystal structure prediction

The basic principle selected for the research to predict the stable nature of structures was to identify the structures with low lattice energy. The crystal phases of temozolomide at global energy minimum are determined from the first principle mentioned in Eq. (1).

The stable crystal structures of temozolomide were identified in the energy landscape plot. The stable hypothetical dense crystal phases of the molecule were generated using MOLPAK software [9]. The algorithm was designed to develop the possible packings by 3-dimensional rotations and repetitions of the parent molecule, incorporating the symmetry operators in all common coordination geometry. The global search was initiated through a gas phase optimization using the Becke 3-Parameter exchange method (B3LYP) [10]. An initial level optimization to minimize the packing density of the hypothetical structures with threshold interactions were selected for PMIN (Packing minimization) optimizations [9] incorporating a repulsion alone potential field called UMD potential [9]. The current research specifically aimed to search the putative crystal forms of temozolomide within the commonly encountered space groups P1, P-1, P2, Pm, Pc, P21, P2/c, P21/m, P2/m, P21/c, Cc, C2, C2/c, Pnn2, Pba2, Pnc2, P221, Pmn21, Pma2, P212121, P21212, Pca21, Pna21, Pnma and Pbca; where 95% [11] of the crystal structures in the database were listed. The PMIN optimization carried out in the global search within the repulsion alone potential field neglected the electron correlation parameters. For the further lattice energy minimization using DMACRYS code [12], the densest hypothetical crystal structures generated from the MOLPAK search were considered, where the structures were subjected to rigid body optimization within a repulsion-dispersion potential field of the form of Eq. (2).where i and k represents the atomic species of 1 and 2 molecules. The FIT potential, calibrated by Williams and Cox [13] was implemented in the lattice energy minimization. The FIT potential was also included with hydrogen atoms attached with nitrogen scaled by Coombes et al. [14]. The GDMA algorithm [15] models electrostatic and intermolecular interactions in the crystal structure of temozolomide from a distributed multipole analysis via MP2/6-31G(d,p) charge density optimization. The thermodynamic stability of the lattice minimized temozolomide molecule was determined from the Ewald summed charge-charge, charge-dipole, dipole-dipole interactions. The second derivative properties of the rigid packings were calculated by the algorithm to justify the mechanical stability of temozolomide molecule. The crystal packings were minimized by removing the symmetrical constraints without achieving the Born criteria resulted in the negative-eigen representations. The unique hypothetical dense crystal packings of temozolomide molecule were ranked to construct the lattice energy landscape.

The lattice energy and density scatterplot (Figure 2) and the energy rank (Table 1) revealed that the predicted crystal structure AM resembles the experimental crystal structure by exhibiting the monoclinic packing of P21/c space group with cell density of 1.587 gm/cm3 at Ulattice = -141.8- kJ/mol. The most favorable crystal packings at 0K along with the experimentally known crystal structure of temozolomide were listed in the lattice energy plot, where comparatively dense region found to be in between -112 kJ/mol to -142 kJ/mol.

Figure 2

Crystal energy landscape showing all possible stable crystal structures of temozolomide along with experimental structure [EXP] in black, Close to real [AM] in blue.

Table 1

List of predicted possible stable crystal structures of temozolomide molecule. The experimental and the predicted, close-to-real polymorphs are in bold.

Crystal	Space	a (Å)	b (Å)	c (Å)	alpha (°)	beta	gamma (°)
ID	group					(°)
AM	P21/c	6.721	7.552	16.101	90	96.1	90
CC	Pbca	7.498	7.752	27.439	90	90	90
CB	Pbca	16.647	7.574	12.646	90	90	90
CA	P-1	6.233	6.119	12.506	119.4	97	75.9
AB	P-1	4.416	9.139	9.9	99.3	91.3	94.9
EXP	P21/n	6.672	7.875	16.02	90	95.3	90
DC	C2/c	19.079	4.563	18.783	90	102.7	90
AH	P21	8.679	6.265	7.661	90	106.9	90
DB	C2	7.715	6.552	16.963	90	105.9	90
FA	P21/c	5.854	25.756	7.855	90	78.4	90
AQ	P212121	14.121	6.554	8.787	90	90	90
CE [Z' = 2]	Pbcn	36.66	7.481	6.241	90	90	90
FC	P21/c	9.796	4.625	17.983	90	80.9	90
AI	P21/c	10.636	4.13	19.359	90	68.7	90
AL	P21/c	7.92	7.011	15.387	90	80.9	90
AY	Pca21	19.399	4.536	8.893	90	90	90
AA	P1	4.232	5.857	8.594	97.1	76.7	81.5
CD[Z' = 2]	P21/c	29.794	7.758	7.101	90	90	90
AZ	P212121	35.714	4.574	4.852	90	90	90
AF	P21	4.819	4.602	17.98	90	92.2	90
AS	Pna21	36.122	4.845	4.572	90	90	90
BF[Z' = 2]	Pm	4.853	36.082	4.569	90	90	90
AP	P21212	4.603	38.241	4.772	90	90	90
DA	Cc	7.316	18.098	9.613	90	139.6	90
BD	Pna21	25.784	6.14	5.109	90	90	90
AV	Pna21	12.988	16.167	4.074	90	90	90
BH	Pca21	18.985	5.057	8.744	90	90	90
EXPb Z’ = 2	P 21/c	17.332	7.351	13.247	90	109.56	90
EXPc Z’ = 2	P -1	8.5	10.004	11.309	99.27	108.86	109.19

Total energy calculation

The total lattice energy Etot of the crystals was approximated with software Crystal Explorer19 [16] on a cluster of molecules surrounding the central one. The N–H and C–H bond distances were constrained to standard values obtained from neutron diffraction studies [17]. The Etot values, were calculated from the sum of scaled electrostatic, polarization, dispersion and repulsion components using the coefficients 1.057, 0.74 and 0.871 respectively using B3LYP/6-31G∗∗ level of theory. For the structures with two molecules per asymmetric unit (Z’ = 2), the energy of the two molecules was averaged. Fingerprint plots of the intermolecular contacts were also generated with program Crystal Explorer19.

Results and discussion

Structure and intermolecular interactions

The temozolomide molecule comprises tetrazine, imidazole and acetamide groups which extends the possibilities of polymorphism. The predicted crystal AM has unit cell parameters close to the real crystal structure and the same space group. The bond lengths, bond angles and dihedral angles of the predicted structure [AM] agree well with that of X-ray geometry [EXP] [5]. There is no significant bond twist for all the bonds of the EXP structure. All the torsion angles in the predicted structures except AM which is close to real one were well agree with the experimentally [EXP] observed structures. The torsion angles of AM structure is varied with that of EXP structure with the minimum and maximum difference of 0.01° and 7.46°. The maximum torsion angle difference is noted for carbonyl (C=O) and methyl (-CH3) attached C3–N2 bond of inner side of tetrazine ring. Figure S1 in supplementary material visualize the torsion angle difference distribution of AM with that of EXP structure. The tetrazine and imidazole rings are almost (with the maximum deviation of 0.02°] in same plane for EXP and this planarity was feebly distracted by 4.7° in AM predicted structure. This planarity was well established from the N6–H6…N3 intra-molecular interaction [Table 2] in both EXP and AM crystal structure. In the solid state confirmation, both the oxygen atoms [O1 and O2] and nitrogen atoms [N4 of tetrazine, N6 of acetamide and N3 of imidazole ring] participated in the intermolecular interactions. The two crystal forms share, in part, similar morphology [Figure 3a, b] and are both stabilized by the dimer interactions framed from N6–H5…O2 strong hydrogen bond [Table 2].

Table 2

Hydrogen bonding details in predicted close to real [AM] and experimental crystal structure [EXP]. Hydrogen bonds common in the two structures are in bold.

Symmetry operation	H-Bond	H…A (Å)	D…A (Å)	D-H … A (˚)
Intra	N6–H6...N3	2.340/2.257	2.756/2.730	103.5/107.1
-x,-y+2,-z/-x-1,-y,-z+2	N6–H5...O2	1.845/1.834	2.836/2.824	166.4/166.0
-x+2,-y+1,-z/-x+1,-y+1,-z+2	C4–H1...O1	2.442/2.370	3.432/3.390	151.3/156.2
x+½,-y+½+1,+z+½/x+½,-y+½,+z-½	C5–H2...N3	2.559/2.700	3.332/3.357	127.6/118.7
-x+½+2,+y+½,-z+½/-x+½+1,+y-½,-z+½+1	C5–H3...O1	2.935/2.935	3.409/3.277	106.7/98.6
-x+½+2,+y+½,-z+½/-x+½+1,+y-½,-z+½+1	C5–H2...O1	2.941/2.713	3.409/3.277	106.4/112.0
-x+½+1,+y-½,-z+½/-x+½,+y+½,-z+½+1	C5–H4...O2	2.564/2.562	3.528/3.571	147.9/154.6
-x+½+1,+y-½,-z+½/-x+½,+y+½,-z+½+1	C5–H4...N4	2.656/2.704	3.464/3.498	130.9/129.8
-/x+1,+y,+z	C5–H3...O2	-/2.838	-/3.532	-/121.9
-/x+1,+y,+z	C5–H3...N4	-/2.760	-/3.571	-/154.6

Figure 3

Overlay of the (a) molecule and (b) supramolecular packing along a axis.

Some hydrogen bonds such as C–H…O, C–H…N and longer distance N–H…N and N–H…O are contributing to the AM and real crystal packings. The intermolecular interactions of the molecule with the neighbouring molecules in the crystal are shown in Figure 4. The Hirshfeld surface analysis gives the details of all the contribution of intermolecular contacts and especially major role played by H...O and H…N contacts. The fingerprint analysis (Figure 5) [18] shows that the N...H/H...N and O...H/H...O contacts give the largest contribution, which are about 28%[EXP]/28.1%[AM] and 25.5%[EXP]/25.9%[AM] respectively. Moreover the H...H contact also renders a significant contribution of 19.7%[EXP]/18.7%[AM] to the crystal packing of the molecule.

Figure 4

View of the strong N–H…O∗ and C–H…O# Hydrogen bond interactions of the molecule in the real crystal. [symmetry operators: ∗-x,-y+2,-z, #-x+2,-y+1,-z].

Figure 5

2D Fingerprint plots of the Hirshfeld surface (a) EXP (b) AM and (c) Difference between EXP and AM.

The fingerprint plot difference between the experimental and predicted molecule was displayed in Figure 5c. In general, the difference were represented by blue and red regions; Blue regions are more intense in the experimental structure and red region highlights for the predicted one. As both the crystal forms prevailed similar interactions with same intense, there is no red region in the difference plot. There is a blue region in the spike at de = 1.1 Å, di = 0.7 Å and de = 0.7 Å, di = 1.1 Å clearly illustrates that there is a greater prevalence of O…H/H…O interactions in experimental structure when compared with that of the predicted one.

From Figure 6, it appears that some of the predicted crystal structures have a two strong hydrogen bonds while the experimental and AM crystal structures have only one. Up to ten predicted crystal structures have two strong hydrogens N–H…O/N involving the two H atoms of the NH2 group. In the real crystal, one NH2 hydrogen atom form one strong N–H…O hydrogen bond and the other one forms two longer distance weak H bonds. The crystal structure FA is characterized by the largest number of hydrogen bonds (2 strong and 7 weak). The BH packing with no strong H-bonds has the smallest energy magnitudes. The number of strong hydrogen bonds shows a correlation coefficient of 0.64 with both -Ulattice and -Eelec values. The crystal packings with two strong H-bonds appear to have generally the largest Eelec lattice electrostatic energy values in magnitude (Figure S4 in supplementary material).

Figure 6

Total number of hydrogen bonds in the different crystal packings. For the Z’ = 2 crystals, the number of H-bonds has been divided by two.

Contact enrichment ratio

The MoPro Viewer [19] package was used to calculate the enrichment ratio which ompares the actual contact with equi-probable theoretical contacts. The interactions; N…C, N…Hc and O… Hc were considered as major interactions as they represent ∼41% of the global Hirshfeld surface and their corresponding enrichment values are listed in Table 3. The strong interaction H5…O2 was represented with E values 1.64 and 1.73 for EXP and AM respectively. The most enriched contact was noted for non-polar hydrogen and oxygen atom and their E values were 1.83 [EXP] and 1.81 [AM]. The calculated values from the contact enrichment ratio reflects similar trend in the interaction energies. It is noted that there is >99% correlation similarity between the intermolecular interaction for predicted crystal structure [AM] and experimental structure [EXP] of temozolomide, which reveals the structures, were equivalent.

Table 3

Contacts types and their enrichment in the predicted [AM] and experimental [EXP] crystal structure of temozolomide.

structure	EXP					AM
atom	Hp	C	N	O	Hc	Hp	C	N	O	Hc
surface %	13.5	23.2	22.9	15.9	24.6	12.9	22.8	23.9	15.9	24.6
Hn	0.3					0.3
C	6.4	5.9			% contacts	6.5	5.7			% contacts
N	8.1	11.1	4.3			8.6	10.9	4.7
O	7.4	8.3	2.1	0.0		7.5	7.9	2.2	0.1
Hc	4.9	7.3	14.9	15.0	4.0	3.4	7.4	15.5	14.9	4.5
Hn	0.16					0.16
C	1.04	1.17			Enrichment	1.11	1.17			Enrichment
N	1.32	1.11	0.86			1.39	1.06	0.87
O	1.64	1.13	0.28	0		1.73	1.1	0.28	0.05
Hc	0.72	0.65	1.33	1.83	0.64	0.51	0.67	1.33	1.81	0.72

Closed and open shell interactions

The electron density and Laplacian of the electron density mappings of the molecules show the covalent bonding regions and lone pairs of the oxygen and nitrogen atoms that are picturized in Figure 7 (a,b). The topological properties of the electron density of the molecule in gas phase lifted from AM is calculated and the same was compared with that from EXP crystal structure in accordance with quantum theory of atoms in molecules (QTAIM) [20]. It is observed that electron density ρ and Laplacian ∇2ρ values show respectively a 93% and 95% correlation coefficient and the values are listed in Table S1 in supplementary material. Among the homogeneous nitrogen-nitrogen bonds, the bond N5–N4 has the highest electron density value ρcp ∼ 3 e/Å3 at the critical point and it is stronger than the N2–N5 bond which has the value of ∼2.2e/Å3. The C5–N2 bond has low ρcp electron density (AM = 1.73e/Å3 and real = 1.66e/Å3). The charges in their bonding region were highly depleted which was shown from less negative Laplacian values (AM = -16.2e/Å5 and real = -12.7e/Å5) in Table S1 (in supplementary material) and found to be weakest among all the C–N bonds with its low ellipticity value (0.04) suggesting the cylindrical nature of the bond.

Figure 7

2D a) Electron density and b) Laplacian of electron density contour plot of temozolomide.

Energy of intermolecular interactions

The stability of the molecules in the experimental and predicted crystal structure was established from C–H…O, N–H…O and C–H….N types of hydrogen bonds. It appears that most other interactions are predominantly stabilized by dispersion energies in predicted crystal structure [AM] and experimental structure (Table 4). Also, the net energies were also calculated from Crystal Explorer 17 using B3LYP/6-31G(d,p) and found the same trend with net electrostatic energies; -98.5 Kj/mol [AM]/-90.1 Kj/mol[EXP] and net dispersion energies; -106.8 kj/mol[AM]/-101.7 Kj/mol[EXP]. From the PIXEL method [21], the intermolecular energy of predicted [AM] and real crystal [EXP] were calculated and are listed in Table 5. Among the interactions, the N6–H5...O2 hydrogen bond is found to be the strongest in the real and predicted AM crystal structures and the corresponding total energy values are -59 kJ/mol and -53 kJ/mol respectively. The similar trend was observed for antipyrine-like derivatives [22]. Figure 8 clearly indicates that the electrostatic, dispersion and total energy frameworks for AM and EXP packed along the a-axis, b-axis and c-axis respectively. This was well reflected from the large and red regions labelled 1 in dnorm maps [Figure 9]. Next to this, the C5–H4...O2 and C5–H4...N4 hydrogen bonds exhibit total energy values of around -33 kJ/mol, and this combined participation was evident from the white regions apart from each other labelled 2 in dnorm maps. The small red regions labelled 3 represents the O1…H1 interactions generating R22(10) motif predicted with total interaction energies of -22.7 kj/mol and -19.7 kj/mol for real [EXP] and predicted [AM] crystal structures respectively.

Table 4

PIXEL lattice energies (kJ/mol) for predicted, AM and real crystal structure, EXP.

term	AM	EXP
Ecoul	-102.3	-97.0
Epol	-30.7	-29.3
EDisp	-111.2	-117.2
ERep	108.3	102.4
Etotal	-135.9	-141.1

Table 5

Interaction energies (kJ/mol) of the molecular pairs, computed with program PIXEL for the real (first lines) and AM unit cell close to real (second lines) polymorphs.

H-Bond	Centroid distance	Ecoul	Epol	EDisp	ERep	Etotal
N6–H5...O2	9.79.6	-72.4-87.5	-24.8-30.2	-19.6-20.3	57.785	-59-53
C4–H1...O1	7.47.5	-20.3-15.5	-5.6-4.6	-13.2-11.8	16.412.1	-22.7-19.7
C5–H2...N3	8.28.4	-4.6-5.2	-2.7-2.2	-12.7-11.5	8.47.4	-11.5-11.6
C5–H3...O1 C5–H2...O1	9.69.5	-3.7-4.0	-1.3-1.3	-6.0-6.2	2.83.6	-8.3-8.0
C5–H4...O2 C5–H4...N4	5.96.0	-20.2-19.1	-5.6-5.5	-19.3-18.8	10.811.4	-34.3-32.1
C5–H3...N4	6.7	-20.9	-6.3	-21.5	16.9	-31.9

Figure 8

Energy frameworks along (a) a-axis, (b) b-axis, (c) c-axis for predicted [AM] and experimental [EXP] structure, showing separate electrostatic (left, red), and dispersion (middle, green) components, and total energy interactions (right, blue). The tube size (scale factor) used in all energy frameworks was 50 and the cut-off was 3.00 kJ/mol.

Figure 9

Hirshfeld surfaces mapped with dnorm for both (a) EXP and (b) AM crystal structures.

The total energy, Etot ​quantified as scaled sum of electrostatic, dispersion, polarization and exchange-repulsion terms from Crystal Explorer program is listed in Table 6 along with Ulattice ​energy. The Etot energy for the real crystal is ordered in third position (−189.4 kJ/mol). The structure AM, close to real, comes, in first position for both Etot = -207.9 kJ/mol and Ulattice = -141.8 kJ/mol.

Table 6

The Ulattice and scaled Etot energy values from DMACRYS and Crystal Explorer software.

Crystal ID	Ulattice Kj/mol	EtotKj/mol
AM	-141.8	-207.9
CC	-140.9	-137.2
CB	-140.7	-158.7
CA	-134.9	-190.3
AB	-132.3	-189.4
EXP	-132.2	-189.4
DC	-130.9	-159.7
AH	-130.9	-135.1
DB	-130.2	-179.2
EXPb Z’ = 2	-129.6	-142.9
FA	-128.0	-178.5
EXPc Z’ = 2	-128.0	-145.5
AQ	-128.0	-127.4
FC	-127.4	-161.8
CE [Z' = 2]	-127.4	-122.9
AI	-127.3	-133.2
AL	-126.7	-137.7
AY	-126.5	-115.2
AA	-125.9	-116.9
CD[Z' = 2]	-125.6	-150.8
AZ	-123.7	-114.0
AF	-123.6	-113.5
AS	-123.1	-111.8
BF[Z' = 2]	-123.0	-91.6
AP	-120.1	-109.3
DA	-119.1	-116.7
BD	-118.4	-109.5
AV	-117.4	-106.1
BH	-112.6	-103.8

Lattice energy of crystal packings

Concerning the Ulattice energy descriptor, the structure AM, close to real, has the largest value (−141.8 kJ/mol) and the real structure is in fourth position.

The CE crystal packing, having two molecules in asymmetric unit, holds the strongest Etot lattice value from CrystalExporer17 at -224.8 kJ/mol. Compared to the real crystal, there are a few additional H-bond interactions prevailing in CE, but the corresponding Ulattice energy value -127.4 k/mol is around ∼10 kJ/mol smaller in magnitude than the real structure.

There is a group of predicted structures which have both weak Ulattice energy (less than -130.9 kJ/mol in magnitude) and weak Etot values (less than -135.1 kJ/mol). They present no strong H-bond (BV, AV) or one (BD, DA, AP, AS, AF, AZ, AY, AA, AQ, AH), like the real structure. . On the whole, the Etot values of the structures correlated with the Ulattice energy values by 47%. The computed Etot energies show variations which are double compared to Ulattice but their average magnitudes are similar (−142 vs. -125 kJ/mol). Globally the Ulattice energies show a good correlation of 0.58 with the lattice Eelec values (Figure S3 in supplementary material) obtained from a transferred multipolar atom model, using UBDB database [23]. The 9 most stable polymorphs according to the Crystalexplorer scaled sum Etot energy are shown in Figure 10. The crystal CE has the strongest Etot value (−238 kJ/mol) and the real crystal appears next within a group formed by (CA, AB, real) crystals for which Etot ≅ 189 kJ/mol. Despite the crystal CE has a strong repulsive energy, this component is strongly attenuated (coefficient 0.618) when the Etot descriptor is considered. Only the EXP polymorph corresponds to one of the crystal structures of TMZ reported in the patent literature [4, 7], which may include some cocrystals. The list of predicted crystal structures in SHELX format were listed in Table S2 in supplementary material.

Figure 10

Scatterplot of lattice energy Ulattice and scaled total energy from CrystalExplorer.

Molecular and packing electrostatic potential

The three-dimensional molecular electrostatic potential (MESP) plot for both experimental and predicted AM structure of temozolomide is shown in Figure 11. The polar nature of the molecule was well established from the positive and negative extrema of the potential on the molecular surface. As expected, in both experimental and predicted AM structure, the strongest negative potential is found for the oxygen atoms O1/O2 and their corresponding potential values were -22.3/-54.4 kcal/mol [EXP] and -22.5/-53.1 kcal/mol [AM] respectively. The maximum positive potential was found near the hydrogen atoms of amine group and their average values were +32.7 kcal/mol [EXP] and +35.0 kcal/mol [AM] respectively. The hydrogen bond donor and acceptor capability was well characterized from potential of oxygen atoms and –NH2 groups.

Figure 11

Front [a] and back [b] view of the molecular surface (electron density isosurface at 0.001 a.u) colored according to the molecular electrostatic potential of the experimental structure. Red coloring: electropositive, blue: electronegative. The positive (VS,max) and negative (VS,min) extrema of potentials are indicated by black and blue spheres, respectively. Values of potential are in kcal/mol/e. The equivalent figure for AM structure is show in Figure S5 in supplementary material.

Figure 12 represents a scatterplot of the inner and outer electrostatic potential V on the Hirshfeld surface for the experimental crystal structure of temozolomide molecule. The points on the surface are distinguished by interaction type. The experimental crystal show a good electrostatic complementarity as a correlation Cvv of -0.81% is seen between the interior and exterior potential. The most electropositive region arise from the strong hydrogen bond interaction [O…Hn and N…Hn] of polar hydrogen atoms, whereas the weaker hydrogen bonds C–H…O and C–H…N leads to less electropositive regions. The most electronegative regions, corresponding to O and N atoms, are H-bond acceptors. The O/N…Hn hydrogen bond acceptor regions appear less electronegative than the O/N…Hc ones as the H-bond distance is smaller: the Hirshfeld surface is closer to the O/N atom nucleus resulting in a less negative potential.

Figure 12

Scatterplot of inner and outer total electrostatic potential V on the Hirshfeld surface for the experimental crystal structure. The contacts related to weak C–H…O and C–H…N hydrogen bonds are in color while those related to a strong hydrogen bond N–H…O or N–H…N are in black.

Contact types in the crystal packings

The difference in contact on the Hirsheld surface between the real crystal and predicted crystals with respect to the proportion of contact types is picturized in Figure 13. In terms of contact surface, the real crystal is more involved in weak hydrogen bonds of N…Hc, O…Hc and C…Hn types, compared to the average of predicted structures.

Figure 13

(a) Average percentage Cxy of actual contacts on the Hirshfeld surface in the predicted crystal packings (in blue). The values for the real packing are shown in black. The errors bars represent the standard deviation in the sample of predicted structures, (b) Percentage of actual contacts in the real and predicted crystal packings, grouped by types.

Strong H-bonds O…Hn and N…Hn are also, on average, more represented in the real crystal. Only the AQ crystal structure has significantly more hydrogen bond surface than the real crystal structure, while the AM close to real structure has slightly less H-bonding. On the other hand, the amount of O…C contacts and hydrophobic interaction Hc…Hc and C…Hc is diminished in the real crystal compared to the average of predictions. Percentage of actual contacts in the experimental and predicted crystal packing are shown in Figure S2 in supplementary material.

Hydrogen propensity analysis

The H-bond propensity [HBP] was calculated for the real and predicted crystal structures to get the possible inter and intramolecular interactions. The propensity calculations were performed from the statistical model established for donor and acceptor pairs, and the training dataset utilized for this model was extracted from Cambridge Structural Database (CSD). The HBP values of possible donor-acceptor pairs of real [EXP] and predicted [AM] structures are listed in Table 7. In general the HBP can take values between 0. and 1.; 0 represents the no likelihood of H-bond formation and 1 indicates the maximal likelihood of H-bond interaction.

Table 7

The H-bond propensity calculation. (∗) represents the observed intermolecular interaction.

Donor	Acceptor	EXP	AM
N6	N3	0.89	0.89
N6	O2	0.81∗	0.81∗
N6	N5	0.80	0.80
N6	N4	0.78	0.78
N6	O1	0.69	0.69

From Table 6, it appears that the atom N6 of acetamide groups was calculated as active donor. On the other hand, the nitrogen atoms [N4 and N5 of tetrazine] and the oxygen atoms [O1 and O2] were calculated as active acceptor atoms in the theoretical calculated H-bond interactions of temozolomide. Among them, the H-bond interactions having highest propensity values [mean value: 0.85] were observed in real [EXP] and predicted is close to real structure [AM], whose mean H-bond propensity value is 0.81.

Figure 14 represents the putative structure landscape for real and predicted temozolomide crystal structures. The crystal structure having the maximum number of H-bond interactions thus having the maximum HBP value were represented by the points situated at the right bottom corner of the landscape plot and the points in the top left represents the unstable polymorphs (Figure S3 in supplementary material].

Figure 14

The putative structure landscape for temozolomide real [EXP] and predicted polymorphs.

Conclusion

Ab initio prediction of the possible stable crystal polymorphs of temozolomide molecule, an oral alkylating pharmaceutical drug, was carried out. By analyzing the lattice energy landscape, unresolved polymorphs with theoretical greater stability were found to exist in the lowest energy regions. The lattice energy landscape was generated via global search through the repulsion alone potential field and incorporated lattice minimization using distributed multipole analysis associated with the hypothetical structures. The polymorph AM shows unit cell parameters close to the experimental crystal form of temozolomide (EXP) and shares some common contacts. The stability of the predicted structure AM was studied by PIXEL energy approach and by evaluating the hydrogen bonds related with the structure. Further studies on the intermolecular interactions were carried out through the Hirshfeld surface analysis. The thermodynamically stable polymorph is not in first position but ranks among the structures with strongest Ulattice and Etot lattice energy descriptors.

Declarations

Author contribution statement

David Stephen Arputharaj: Conceived and designed the experiments; Performed the experiments; Wrote the paper.

Meenashi Rajasekaran, Christian Jelsch: Analyzed and interpreted the data.

Saravanan Kandasamy, Abdullah G Al-Sehemi: Contributed reagents, materials, analysis tools or data.

Funding statement

TThis research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Data availability statement

Data will be made available on request.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.