Monoclinic Paracetamol vs . Paracetamol-4 , 4 ′-Bipyridine Co-Crystal ; What Is the Difference ? A Charge Density Study

Paracetamol (PCM) has two well-documented polymorphic forms at room temperature; monoclinic Form I is more stable than the other orthorhombic Form II. Form II exhibits improved tabletting properties compared to Form I due to low shearing forces; however, difficulties in its manufacture have limited its use in industrial manufacture. Previous studies have found that the introduction of a co-former to form co-crystals would allow the PCM molecule to exist in a conformation similar to that of the orthorhombic form while being more stable at room temperature. Experimental charge density analysis of the paracetamol-4,4′-bipyridine (PCM-44BP) co-crystal system, and its constituent molecules, has been carried out to examine the forces that drive the formation and stabilisation of the co-crystal, while allowing PCM to maintain a packing motif similar to that found in Form II. It is hoped studies on this well-known compound will help apply the knowledge gained to other drug molecules that are less successful. The PCM molecules in the co-crystal were found to exhibit similar packing motifs to that found in Form I, however, intercalation of the 44BP molecule between the PCM layers resulted in a shallower angle between molecular planes, which could result in the required lateral shear. Topological analysis identified more weak interactions in the co-crystal compared to the individual molecules, thus allowing for greater stability as evidenced by the lattice energies. Weak interactions in the PCM-44BP co-crystal were found to range in strength from 4.08–84.33 kJ·mol−1, and this variety allowed the PCM-44BP planes to be held together, while a weak π–π interaction (15.14 kJ·mol−1) allowed lateral shear to occur, thus mimicking the planes found in Form II PCM and offering the possibility of improved tabletting properties. A comparison of integrated atomic charges between partitions of the PCM molecules in the single and co-crystal found that the hydroxyl and amide groups were involved in greater hydrogen bonding in the co-crystal, resulting in a charge redistribution across the molecule evidenced by a larger molecular dipole moment (μ = 12.34D). These findings, in addition to the co-crystal having the largest lattice energy, form a potential basis with which to predict that the co-crystal exhibits improved solubility and stability profiles. It is anticipated that these findings will contribute to improvements in the formulation and other physical properties of PCM and other pharmaceutical compounds.


Introduction
Crystal engineering refers to the rational design of solids by inducing the reproducible formation of weak interactions between pairs of functional groups [1]. This field has recently undergone a resurgence because of increased interest from pharmaceutical companies and research institutions. This coincides with a recent reduction in the number of novel drugs approved by the FDA and other global regulatory agencies relative to the early 1990s [2].
Due to the excessive costs associated with bringing a new drug from the design stage through to the market, many pharmaceutical companies have opted to revisit older drug candidates that did not make it to market and to seek to improve their physicochemical properties. This approach encompasses a very large number of drugs as over 90% of drug candidates fail due to poor physicochemical properties, resulting in sub-optimal in vivo performance or a requirement for the development of inconvenient dosage forms [3]. Incorporating these drugs into co-crystal systems has emerged as an innovative and effective means of improving the physicochemical properties of these agents without any requirements for covalent modifications [4]. This is especially relevant in cases where salt formation is not a viable option (i.e., in cases where there is no ionisable group on the drug). There continues to be much debate regarding the definition of co-crystals, including whether solvates and hydrates are included in the same definition: the FDA has defined co-crystals as "solids that are crystalline materials composed of two or more molecules in the same crystal lattice" with the interactions between them governed by "non-ionic interactions" [5].
The use of co-crystals to modulate the physicochemical properties of an active pharmaceutical ingredient (API) has seen some success in the form of dapagliflozin, the sodium/glucose cotransporter 2 (SGLT2) inhibitor for the treatment of type 2 diabetes (T2D) [6]. However, the approach is still in its infancy with current approaches giving unpredictable results. Paracetamol (PCM) has been extensively studied both individually and as part of co-crystal systems; it is known that PCM exists in three forms with the monoclinic Form I being the most thermodynamically stable. Two other, orthorhombic, metastable forms require more sophisticated growth conditions including melting and anhydrous nitrogen conditions for Form III [7][8][9]. Although PCM is well known, the stable monoclinic Form I is known to have very poor compressibility and tabletability properties, while the orthorhombic metastable Form II exhibits improved tabletability [10,11]. Commercial production of PCM tablets currently utilizes Form I PCM in combination with gelatine or starch derivatives to improve its tableting properties [7]. The commercial success of PCM Form I has rendered the poor tabletability a non-issue. However, it is hoped studies on such a well-known drug molecule, resulting in increased understanding of the relationship between different polymorphic forms and their physical properties may lay the fundamentals for applying this knowledge to less successful drugs. Different packing motifs in the PCM polymorphs are believed to be one of the reasons why they exhibit different physical properties with Form I forming zig-zag shapes via O-H···O=C hydrogen bonds and Form II being organised into parallel sheets and slip planes, which make it easier to form into tablets via N-H···O and O-H···O intermolecular hydrogen bonds [12]. It has been established that Form II can only be produced reliably from the melt, and many research groups have thus tried to co-crystallise PCM in attempts to replicate the sheet packing motif and slip planes present in Form II. This has been achieved with limited success; a variety of co-formers have been used, with only three forming the desired motif [12,13]. As epitomised by this example alone, the brute force approach to developing co-crystals with the desired physical properties is no longer feasible, and a more targeted approach is required. Greater insights into the driving forces behind the formation of weak interactions, such as hydrogen bonds, van der Waals and π-π interactions, are required to develop a more streamlined and information-driven approach to developing pharmaceutically-desirable co-crystals, as well as to help develop the field of crystal engineering.
High resolution X-ray diffraction can be used to determine the electron density distribution (EDD) in crystalline systems. Here, we present a comparative analysis of the EDD between PCM Form I (1), 4,4 -bipyridine (44BP) (2)

Methods
Raw materials of (1) and (2) were purchased from Sigma Aldrich and used without further purification. Crystals for (1) were grown via slow evaporation from ethanol. Crystals for (2) were as received from Sigma Aldrich. The co-crystal (3) was grown by dissolving equimolar amounts of (1) and (2) individually, in ethanol, before mixing, and the solution was allowed to evaporate.

Computational Methods
Gas phase, single point (SP) calculations were performed on (1), (2) and (3) with geometry taken from the high-order experimental coordinates. All theoretical calculations were performed with the Gaussian 09 suite [14] at the 6-311++G** level of theory for all structures. All calculations used the three-parameter hybrid exchange function developed by Becke [15] in conjunction (vide supra) with the exchange correlation potential, corrected via gradient developed by Lee et al. [16] together with the long-range correction proposed by Tawada et al. [17,18] (CAM-B3LYP). Analysis of the topology of electron density from the experimental model was performed via the properties module of XD, XDPROP [19], while analysis of the electron density for the theoretical densities was performed using the AIMALL [20] package.
Details on the collection, integration and reduction of data can be found in the Supporting Information. The multipole refinement procedure used has been reported in previous publications [21][22][23]. Refer to Table 1 for selected crystallographic information from the independent atom model (IAM) and multipole (EXP) refinements.

Methods
Raw materials of (1) and (2) were purchased from Sigma Aldrich and used without further purification. Crystals for (1) were grown via slow evaporation from ethanol. Crystals for (2) were as received from Sigma Aldrich. The co-crystal (3) was grown by dissolving equimolar amounts of (1) and (2) individually, in ethanol, before mixing, and the solution was allowed to evaporate.

Computational Methods
Gas phase, single point (SP) calculations were performed on (1), (2) and (3) with geometry taken from the high-order experimental coordinates. All theoretical calculations were performed with the Gaussian 09 suite [14] at the 6-311++G** level of theory for all structures. All calculations used the three-parameter hybrid exchange function developed by Becke [15] in conjunction (vide supra) with the exchange correlation potential, corrected via gradient developed by Lee et al. [16] together with the long-range correction proposed by Tawada et al. [17,18] (CAM-B3LYP). Analysis of the topology of electron density from the experimental model was performed via the properties module of XD, XDPROP [19], while analysis of the electron density for the theoretical densities was performed using the AIMALL [20] package.
Details on the collection, integration and reduction of data can be found in the Supporting Information. The multipole refinement procedure used has been reported in previous publications [21][22][23]. Refer to Table 1 for selected crystallographic information from the independent atom model (IAM) and multipole (EXP) refinements.

Geometry
Bond lengths and angles for all experimental structures were obtained from the multipole model (MM) refinement output. For all of the structures, the X-ray structure was in excellent agreement with results reported in the literature with mean differences in bond lengths of 0.002, 0.044 and 0.001 Å for (1), (2) and (3), respectively. Similarly, the mean differences in angles were 0.086, 0.039 and 0.043 • . Refer to Tables S4-S21 (ESI) for a full list of bond lengths and angles. The crystal structure of (1) contained one molecule of PCM in the asymmetric unit. The asymmetric units of (2) and (3) both contained more than one molecule with (2) containing two 44BP molecules in a T-shaped arrangement and (3) containing two molecules each of PCM and 44BP in an alternating pattern. Refer to  for ORTEP diagrams and labelling schemes.     (2). Thermal ellipsoids are shown at the 50% probability level [24]. Figure 3. ORTEP diagram of 44BP (2). Thermal ellipsoids are shown at the 50% probability level [24].  . Thermal ellipsoids are shown at the 50% probability level [24].
The conformational geometry of PCM was found to vary between the individual and co-crystal forms. The C(05)-C(04)-N(01)-C(07) torsion angle highlighting the orientation of the acetamide group relative to the aromatic ring in (1) was found to be 22.69°. The analogous angles in (3) were found to be 10.28 and 14.49° for the two crystallographically independent PCM molecules. Figure 5 shows the packing of PCM in both (1) and (3). Figure 6 shows the differing conformation of the PCM molecules in (1) and (3) due to the differing torsion angles. The zig-zag packing shape previously described by Ahmed et al. [25] in the monoclinic form of PCM is clearly present in both (1) and (3) and explains why (3) exhibits worse tabletting properties than other PCM co-crystals [13]. However, further analysis of the angle between PCM planes reveals a much larger angle between PCM packing planes (116.91 vs. 149.81° for (1) in the C(04)-C(07)-C(02) plane and (3) in the C(07′)-C(08′)-N(01) plane).
These differences in PCM packing can be attributed to the conformation of the 44BP molecules . Thermal ellipsoids are shown at the 50% probability level [24].
The conformational geometry of PCM was found to vary between the individual and co-crystal forms. The C(05)-C(04)-N(01)-C(07) torsion angle highlighting the orientation of the acetamide group relative to the aromatic ring in (1) was found to be 22.69 • . The analogous angles in (3) were found to be 10.28 and 14.49 • for the two crystallographically independent PCM molecules. Figure 5 shows the packing of PCM in both (1) and (3). Figure 6 shows the differing conformation of the PCM molecules in (1) and (3) due to the differing torsion angles. The zig-zag packing shape previously described by Ahmed et al. [25] in the monoclinic form of PCM is clearly present in both (1) and (3) and explains why (3) exhibits worse tabletting properties than other PCM co-crystals [13]. However, further analysis of the angle between PCM planes reveals a much larger angle between PCM packing planes (116.91 vs. 149.81 • for (1) in the C(04)-C(07)-C(02) plane and (3) in the C(07 )-C(08 )-N(01) plane).
These differences in PCM packing can be attributed to the conformation of the 44BP molecules in (3). There are no major differences between the conformations of the 44BP molecules in (2) and (3) with the dihedral angles (C(4)-C(5)-C(6)-C(7)) being 16.6 and −34.6 • in (2) and 32.4 and 19.9 • in (3). Previous studies have shown that the minimum conformation energy structure of 44BP corresponds to a dihedral angle of 38.6 • [26]. The conformation of the second 44BP molecule in both (1) and (3) resulting in a higher energy state is in accordance with previous observations that 44BP molecules are highly influenced by surrounding interactions such as crystal packing effects [27]. Based on the previous discussion, it can be said that the intercalation of the 44BP molecules into the PCM packing structure in (3) causes the increased angle between PCM packing planes when compared to (1) and results in improved tableting performance over the single constituent crystal.

Topological Analysis
Topological analysis based on Bader's Atoms in Molecules (AIM) [30] theory were carried out on both the EXP and SP models. Completeness of analysis was ensured by ensuring satisfaction of

Topological Analysis
Topological analysis based on Bader's Atoms in Molecules (AIM) [30] theory were carried out on both the EXP and SP models. Completeness of analysis was ensured by ensuring satisfaction of

Topological Analysis
Topological analysis based on Bader's Atoms in Molecules (AIM) [30] theory were carried out on both the EXP and SP models. Completeness of analysis was ensured by ensuring satisfaction of the Poincaré-Hopf or its crystalline equivalent, the Morse relationship [31]. There was very little difference in the topological parameters obtained from the EXP and SP models as seen through the small differences in ρ bcp and ∇ 2 ρ bcp. Refer to Tables S22-S27 (ESI) for a full list of BCP and ring critical points (RCPs) found in our analysis.
For (1), mean differences of 0.044 eÅ −3 and 4.26 eÅ −5 were found for ρ bcp and ∇ 2 ρ bcp , respectively. For (2) and (3), mean differences of 0.053 and 0.020 e Å −3 and 1.68 and 3.76 eÅ −5 were found for ρ bcp and ∇ 2 ρ bcp , respectively. For (1), the largest discrepancies between EXP and SP were in ∇ 2 ρ bcp for the O(01)-H(01), N(01)-H(01A) and O(02)-C(07) bonds. For the first two bonds, the EXP ∇ 2 ρ bcp was greater than the SP value, while the converse was true for the third bond. In all cases, both models report open shell interactions as expected of covalent bonds; however, the differences in magnitude can be attributed to the first two bonds being hydrogen bond donors, while the third is a hydrogen bond acceptor. Similar trends can also be seen in the analogous bonds in (3), which reiterates the findings discussed above that the 44BP molecules seem to intercalate between the PCM molecules and do not interfere with the intermolecular bonding between those molecules. There were no significant differences between the EXP and SP models for (2).

Hydrogen Bonds
There were no clear trends present in the hydrogen bonds found for all systems from a geometrical perspective. One and two intramolecular bonds were found in (1) and (3), respectively. In each case, the bond was located in the PCM molecule between the carbonyl oxygen and aromatic hydrogen (C(05)-H(05)···O(02)). The hydrogen bond acceptor atoms involved in intramolecular bonding in (2) were also involved in similar hydrogen bonding arrangements in (3).
The four hydrogen bonds in (1) all have E HB > 20 kJ·mol −1 . The intramolecular bond is the weakest of the four, while the remaining bonds are all above 50 kJ·mol −1 . The three intermolecular hydrogen bonds spread out in opposite directions to stabilise the two PCM molecules perpendicular to the original. An approximate binding energy of 190 kJ·mol −1 in this area resulting in the formation of a perpendicular plane of PCM is responsible for the zig-zag packing motif and poor compressibility of the monoclinic form of PCM, confirming the observations reported previously [13]. Figure 7 shows the Laplacian of the intermolecular hydrogen bonds, which bond the PCM molecules in (1) together into the zig-zag packing shape. The crystal structure of (2) is held together by weak interactions with all hydrogen bonds being less than 20 kJ mol −1 in strength. The hydrogen bonds between N(2′) and H(10′) and H(1′) help the structure to maintain the T-shaped arrangement in the asymmetric unit. The two bonds however are not of equal strength, with the N(2)···H(10′) being stronger than the N(2)···H(1′) bond (18.8 kJ mol −1 vs. 8.36 kJ mol −1 ). Previous studies [35,36] have examined the concept of lone pair (LP) position as a predictor of hydrogen bond directionality and have found that the closer the value of the DHA (donor-hydrogen-acceptor) and DH-LP (donor-hydrogen-lone pair) angles, the more likely a hydrogen bond is to form. Furthermore, the DH-LP angle may also be reflective of the direction of polarisation of the valence shell of the acceptor atom and may have a direct effect on bond strength. The DH-LP angles for the C(10′)-H(10′)···N(2) and C(1′)-H(1′)···N(2) bonds are 150.8 and 143.4° degrees, respectively. This is also seen in Figure 8a,c where the nitrogen lone pair is in a more linear alignment towards H(10′). The difference in DH-LP angle may be the cause of the disparity in bond strength. The difference in strength between these two stabilising interactions may contribute to the presence of two 44BP molecules in the asymmetric unit with differing geometries.
Although the intermolecular bonds present in (2) were not conserved in (3), the nitrogen atoms play a similar role and form very similar bonding motifs. This can be seen in Figure 8b,e, where analogous nitrogen atoms (N(2) in (2) and N(2′) in (3) both form bifurcated hydrogen bonds and can be seen both EXP and SP models. Figure 8e shows a third hydrogen bond between N(2′) and H(03′) in the SP model for (3). This bond was not found in the topological analysis of the EXP model. Interestingly, the DH-LP angles for the N(2′)···H(03′) and N(2′)···H(08E) bonds are very similar with values of 128.6 and 132.8°, respectively. Visual analysis of the polarisation of the N(2′) atom in Figure 8d shows that for the EXP model, the lone pair is almost directly facing H(01′) with a slight deviation towards H(08E), whereas in Figure 8e, the polarization of the same atom appears to have some spread towards H(03′). This difference in polarisation is also reflected in the difference in bond strengths. In the EXP model, both N(2′)···H(08E) and N(2′)···H (01)   The crystal structure of (2) is held together by weak interactions with all hydrogen bonds being less than 20 kJ·mol −1 in strength. The hydrogen bonds between N(2 ) and H(10 ) and H(1 ) help the structure to maintain the T-shaped arrangement in the asymmetric unit. The two bonds however are not of equal strength, with the N(2)···H(10 ) being stronger than the N(2)···H(1 ) bond (18.8 kJ·mol −1 vs. 8.36 kJ·mol −1 ). Previous studies [35,36] have examined the concept of lone pair (LP) position as a predictor of hydrogen bond directionality and have found that the closer the value of the DHA (donor-hydrogen-acceptor) and DH-LP (donor-hydrogen-lone pair) angles, the more likely a hydrogen bond is to form. Furthermore, the DH-LP angle may also be reflective of the direction of polarisation of the valence shell of the acceptor atom and may have a direct effect on bond strength. The DH-LP angles for the C(10 )-H(10 )···N(2) and C(1 )-H(1 )···N(2) bonds are 150.8 and 143.4 • degrees, respectively. This is also seen in Figure 8a,c where the nitrogen lone pair is in a more linear alignment towards H(10 ). The difference in DH-LP angle may be the cause of the disparity in bond strength. The difference in strength between these two stabilising interactions may contribute to the presence of two 44BP molecules in the asymmetric unit with differing geometries.
Although the intermolecular bonds present in (2) were not conserved in (3), the nitrogen atoms play a similar role and form very similar bonding motifs. This can be seen in Figure 8b,e, where analogous nitrogen atoms (N(2) in (2) and N(2 ) in (3) both form bifurcated hydrogen bonds and can be seen both EXP and SP models.  relatively planar conformation, with the slightly stronger bonds in (3) compared to (1) due to the larger torsion angle between the aromatic ring and acetamide group as discussed above. Examination of the packing motifs in (3) show that the 44BP molecules stack on top of each other and form a bridge between two PCM molecules as seen in Figure 9. The O(01)-H(01A)···O(02 ) bond (84.33 kJ·mol −1 ) is the strongest hydrogen bond in (3) and links the different packing motifs together by joining the hydroxyl group of one PCM molecule to the acetamide group of another. One nitrogen on each of the 44BP molecules (N(2) and N(2 ) forms moderate strength bonds (approximately 45 and 60 kJ·mol −1 , respectively), while the remaining nitrogen forms weak bonds. N(2 ) forms two hydrogen bonds to its closest PCM neighbour. The two bonds formed by N(2 ) along with the bond formed by N(1 ) result in a twist in the C(5 )-C(6 ) axis to a dihedral angle of 19.88 • , while previous literature has found that the minimum energy state of 44BP corresponds to a dihedral angle of 38.6 • [26]. The second 44BP molecule, which is not as involved in hydrogen bonding, has a conformation closer to the energy minimum with a dihedral angle of 32.4 • . The π-π interactions are present between the stacked 44BP molecules shown in Figure 9, and the interaction strength was found to be 15.14 kJ·mol −1 using the method developed by Waller et al. [37] on the SP model. This interaction can be considered weak. Overall, the hydrogen bonds discussed above ensure that the planes of PCM and 44BP molecules are held together to maintain physical stability; however, the weak interaction between 44BP molecules allows lateral shear to occur, potentially resulting in improved compressibility and tabletability of the co-crystal like the properties reported for form II PCM.
An inverse relationship is known to exist between bond strength and length [38]. This relationship was also investigated, specifically the distances between the hydrogen and acceptor atoms of the bond and the bond critical point (bcp) (d H···bcp and d A···bcp ). The relationship was found to be true for the hydrogen bonds found across all systems.  [26]. The second 44BP molecule, which is not as involved in hydrogen bonding, has a conformation closer to the energy minimum with a dihedral angle of 32.4°. The π-π interactions are present between the stacked 44BP molecules shown in Figure 9, and the interaction strength was found to be 15.14 kJ mol −1 using the method developed by Waller et al. [37] on the SP model. This interaction can be considered weak. Overall, the hydrogen bonds discussed above ensure that the planes of PCM and 44BP molecules are held together to maintain physical stability; however, the weak interaction between 44BP molecules allows lateral shear to occur, potentially resulting in improved compressibility and tabletability of the co-crystal like the properties reported for form II PCM. An inverse relationship is known to exist between bond strength and length [38]. This relationship was also investigated, specifically the distances between the hydrogen and acceptor atoms of the bond and the bond critical point (bcp) (dH···bcp and dA···bcp). The relationship was found to be true for the hydrogen bonds found across all systems.

Hirshfeld Surfaces Analysis
Three-dimensional Hirshfeld surfaces and corresponding two-dimensional fingerprint plots were generated using the CrystalExplorer program [39]. Hirshfeld surfaces provide a visual representation of the electron density around a molecule, illustrating the extent to which the electrons in the space surrounding the molecules can be attributed to the molecule. Hirshfeld surfaces for (1), (2) Figures S4-S6, while Figure 10 shows selected fingerprint plots for (3). These plots are an extension of the Hirshfeld surface and provide information on the distance of the interacting nuclei from the surface. Values of de and di are defined as the distance from the Hirshfeld surface to the nearest external and internal nucleus, respectively. The parameter dnorm is defined as the normalised contact distance with reference to de, di and van der Waals radii of the atoms [40]. Unsurprisingly, analysis of the fingerprint plots for all systems shows that H···H and C···C bonds contribute the most towards the Hirshfeld surface in the form of van der Waals interactions. The remainder of the fingerprint plots are composed of the weak interactions in the system. The O···H and N···H interactions also play a role in maintaining the structure of (1) and (2), accounting for 26.2% and 20.7% of the fingerprint plot, respectively. Both of these interactions contribute evenly to the fingerprint plot of (3), with 11.9% and 9.5%, respectively. Similar plots for (1) and (2) can be found in the Supplementary Information in Figures S7 and S8, respectively.

Hirshfeld Surfaces Analysis
Three-dimensional Hirshfeld surfaces and corresponding two-dimensional fingerprint plots were generated using the CrystalExplorer program [39]. Hirshfeld surfaces provide a visual representation of the electron density around a molecule, illustrating the extent to which the electrons in the space surrounding the molecules can be attributed to the molecule. Hirshfeld surfaces for (1), (2) and (3) can be found in the Supplementary Information in Figures S4-S6, while Figure 10 shows selected fingerprint plots for (3). These plots are an extension of the Hirshfeld surface and provide information on the distance of the interacting nuclei from the surface. Values of d e and d i are defined as the distance from the Hirshfeld surface to the nearest external and internal nucleus, respectively. The parameter d norm is defined as the normalised contact distance with reference to d e , d i and van der Waals radii of the atoms [40]. Unsurprisingly, analysis of the fingerprint plots for all systems shows that H···H and C···C bonds contribute the most towards the Hirshfeld surface in the form of van der Waals interactions. The remainder of the fingerprint plots are composed of the weak interactions in the system. The O···H and N···H interactions also play a role in maintaining the structure of (1) and (2), accounting for 26.2% and 20.7% of the fingerprint plot, respectively. Both of these interactions contribute evenly to the fingerprint plot of (3), with 11.9% and 9.5%, respectively. Similar plots for (1) and (2) can be found in the Supplementary Information in Figures S7 and S8, respectively.

Atomic Charges
Integrated atomic basin charges were also determined from topological analysis of (1), (2) and (3). There was excellent agreement between the charges obtained from the EXP and SP models. Across all systems, the largest differences were seen in the nitrogen atoms in the 44BP molecules in (2) and (3). The EXP model consistently underestimates the charge on these atoms by approximately 0.3e. This can be attributed to their involvement in intermolecular hydrogen bonding often as acceptors as previously discussed. The gas phase nature of the SP calculations along with only the asymmetric unit being evaluated removing supramolecular effects such as crystal packing may also affect the charge distribution. This is especially relevant in the case of the 44BP molecule, as it is well known that the conformation is heavily influenced by its supramolecular environment [27,41]. Similar analysis of the atomic charges between atoms in (1) and (2) and their analogous atoms in (3) was insignificant, and no differences of note were found.
The PCM and 44BP molecules in all structures were partitioned into distinct chemical groups to examine the changes in charge across each part of the molecule. PCM was separated into three sections; the -OH, aromatic ring and -CONH2CH3 groups. The -CONH2CH3 groups underwent the largest change going from having a net charge of −0.39e in (1) to −0.13 and −0.17e in (3). This is due to the group being significantly more involved in the hydrogen bond in the co-crystal where it contributes as a hydrogen bond acceptor four times and five times as a donor. The hydroxyl group also follows this trend, being involved in four hydrogen bonds including the O(01)-H(01A)···O(02′) bond, which was the strongest in (3) with a strength of 84.33 kJ mol −1 . As required, the aromatic ring

Atomic Charges
Integrated atomic basin charges were also determined from topological analysis of (1), (2) and (3). There was excellent agreement between the charges obtained from the EXP and SP models. Across all systems, the largest differences were seen in the nitrogen atoms in the 44BP molecules in (2) and (3). The EXP model consistently underestimates the charge on these atoms by approximately 0.3e. This can be attributed to their involvement in intermolecular hydrogen bonding often as acceptors as previously discussed. The gas phase nature of the SP calculations along with only the asymmetric unit being evaluated removing supramolecular effects such as crystal packing may also affect the charge distribution. This is especially relevant in the case of the 44BP molecule, as it is well known that the conformation is heavily influenced by its supramolecular environment [27,41]. Similar analysis of the atomic charges between atoms in (1) and (2) and their analogous atoms in (3) was insignificant, and no differences of note were found.
The PCM and 44BP molecules in all structures were partitioned into distinct chemical groups to examine the changes in charge across each part of the molecule. PCM was separated into three sections; the -OH, aromatic ring and -CONH 2 CH 3 groups. The -CONH 2 CH 3 groups underwent the largest change going from having a net charge of −0.39e in (1) to −0.13 and −0.17e in (3). This is due to the group being significantly more involved in the hydrogen bond in the co-crystal where it contributes as a hydrogen bond acceptor four times and five times as a donor. The hydroxyl group also follows this trend, being involved in four hydrogen bonds including the O(01)-H(01A)···O(02 ) bond, which was the strongest in (3) with a strength of 84.33 kJ·mol −1 . As required, the aromatic ring in PCM became more positive to compensate. The 44BP molecule was separated along its axis of symmetry. Due to the 44BP molecules being involved in similar bonding motifs in (2) and (3), it is unsurprising that there was very little difference in charge with the largest difference being 0.11e. Details of each moiety and their respective charges can be found in Table 3. The complete list of atomic charges can be found in Tables S33-S35 in the Supplementary Information. in PCM became more positive to compensate. The 44BP molecule was separated along its axis of symmetry. Due to the 44BP molecules being involved in similar bonding motifs in (2) and (3), it is unsurprising that there was very little difference in charge with the largest difference being 0.11e. Details of each moiety and their respective charges can be found in Table 3. The complete list of atomic charges can be found in Tables S33-S35 in the Supplementary Information.

Electrostatic Potential
The molecular electrostatic potential (MEP) is an integral part in linking the results found in EDD studies such as this to real physical outcomes such as atom reactivities and how electrons move when weak interactions are formed. The MEP maps shown in Figures 11-13 are a visualisation of these properties, and a comparison of these maps for the molecules in their individual and co-crystalline states provides insights into which groups are involved in the formation of weak interactions and the impetus behind them to form. The MEP maps have been mapped onto an isosurface of ρ, which have been plotted on the same scale for comparability. Visual inspection of the maps for (1), (2) and (3) did not yield any unusual features. The most electronegative regions are located near the hetero atoms, and the most electropositive regions were found near the hydrogen atoms. This is expected as most weak interactions found were traditional hydrogen bonds, and these areas correspond to regions where much of the intermolecular bonding occurs in the form of hydrogen bonds. Further analysis of the map of (3) shows excellent complementarity between bond donors and acceptors with them being more electropositive and electronegative, respectively. It was also found that larger differences in electrostatic potential between bond donor/acceptor pairs corresponded to stronger hydrogen bonds, as discussed above. Additionally, the atoms on the periphery of the asymmetric unit with notable electrostatic potential charges were the same atoms involved in intermolecular bonding as discussed previously. The findings from the analysis of the MEP maps in Figures 11-13 validate the future use of electrostatic potential as a predictor of hydrogen bonds and other weak interactions as opposed to the current geometrical definitions used today. in PCM became more positive to compensate. The 44BP molecule was separated along its axis of symmetry. Due to the 44BP molecules being involved in similar bonding motifs in (2) and (3), it is unsurprising that there was very little difference in charge with the largest difference being 0.11e. Details of each moiety and their respective charges can be found in Table 3. The complete list of atomic charges can be found in Tables S33-S35 in the Supplementary Information. Table 3. Atomic charges (e) from multipole refinement.

Electrostatic Potential
The molecular electrostatic potential (MEP) is an integral part in linking the results found in EDD studies such as this to real physical outcomes such as atom reactivities and how electrons move when weak interactions are formed. The MEP maps shown in Figures 11-13 are a visualisation of these properties, and a comparison of these maps for the molecules in their individual and co-crystalline states provides insights into which groups are involved in the formation of weak interactions and the impetus behind them to form. The MEP maps have been mapped onto an isosurface of ρ, which have been plotted on the same scale for comparability. Visual inspection of the maps for (1), (2) and (3) did not yield any unusual features. The most electronegative regions are located near the hetero atoms, and the most electropositive regions were found near the hydrogen atoms. This is expected as most weak interactions found were traditional hydrogen bonds, and these areas correspond to regions where much of the intermolecular bonding occurs in the form of hydrogen bonds. Further analysis of the map of (3) shows excellent complementarity between bond donors and acceptors with them being more electropositive and electronegative, respectively. It was also found that larger differences in electrostatic potential between bond donor/acceptor pairs corresponded to stronger hydrogen bonds, as discussed above. Additionally, the atoms on the periphery of the asymmetric unit with notable electrostatic potential charges were the same atoms involved in intermolecular bonding as discussed previously. The findings from the analysis of the MEP maps in Figures 11-13 validate the future use of electrostatic potential as a predictor of hydrogen bonds and other weak interactions as opposed to the current geometrical definitions used today. in PCM became more positive to compensate. The 44BP molecule was separated along its axis of symmetry. Due to the 44BP molecules being involved in similar bonding motifs in (2) and (3), it is unsurprising that there was very little difference in charge with the largest difference being 0.11e. Details of each moiety and their respective charges can be found in Table 3. The complete list of atomic charges can be found in Tables S33-S35 in the Supplementary Information. Table 3. Atomic charges (e) from multipole refinement.

Electrostatic Potential
The molecular electrostatic potential (MEP) is an integral part in linking the results found in EDD studies such as this to real physical outcomes such as atom reactivities and how electrons move when weak interactions are formed. The MEP maps shown in Figures 11-13 are a visualisation of these properties, and a comparison of these maps for the molecules in their individual and co-crystalline states provides insights into which groups are involved in the formation of weak interactions and the impetus behind them to form. The MEP maps have been mapped onto an isosurface of ρ, which have been plotted on the same scale for comparability. Visual inspection of the maps for (1), (2) and (3) did not yield any unusual features. The most electronegative regions are located near the hetero atoms, and the most electropositive regions were found near the hydrogen atoms. This is expected as most weak interactions found were traditional hydrogen bonds, and these areas correspond to regions where much of the intermolecular bonding occurs in the form of hydrogen bonds. Further analysis of the map of (3) shows excellent complementarity between bond donors and acceptors with them being more electropositive and electronegative, respectively. It was also found that larger differences in electrostatic potential between bond donor/acceptor pairs corresponded to stronger hydrogen bonds, as discussed above. Additionally, the atoms on the periphery of the asymmetric unit with notable electrostatic potential charges were the same atoms involved in intermolecular bonding as discussed previously. The findings from the analysis of the MEP maps in Figures 11-13 validate the future use of electrostatic potential as a predictor of hydrogen bonds and other weak interactions as opposed to the current geometrical definitions used today. in PCM became more positive to compensate. The 44BP molecule was separated along its axis of symmetry. Due to the 44BP molecules being involved in similar bonding motifs in (2) and (3), it is unsurprising that there was very little difference in charge with the largest difference being 0.11e. Details of each moiety and their respective charges can be found in Table 3. The complete list of atomic charges can be found in Tables S33-S35 in the Supplementary Information.

Electrostatic Potential
The molecular electrostatic potential (MEP) is an integral part in linking the results found in EDD studies such as this to real physical outcomes such as atom reactivities and how electrons move when weak interactions are formed. The MEP maps shown in Figures 11-13 are a visualisation of these properties, and a comparison of these maps for the molecules in their individual and co-crystalline states provides insights into which groups are involved in the formation of weak interactions and the impetus behind them to form. The MEP maps have been mapped onto an isosurface of ρ, which have been plotted on the same scale for comparability. Visual inspection of the maps for (1), (2) and (3) did not yield any unusual features. The most electronegative regions are located near the hetero atoms, and the most electropositive regions were found near the hydrogen atoms. This is expected as most weak interactions found were traditional hydrogen bonds, and these areas correspond to regions where much of the intermolecular bonding occurs in the form of hydrogen bonds. Further analysis of the map of (3) shows excellent complementarity between bond donors and acceptors with them being more electropositive and electronegative, respectively. It was also found that larger differences in electrostatic potential between bond donor/acceptor pairs corresponded to stronger hydrogen bonds, as discussed above. Additionally, the atoms on the periphery of the asymmetric unit with notable electrostatic potential charges were the same atoms involved in intermolecular bonding as discussed previously. The findings from the analysis of the MEP maps in Figures 11-13 validate the future use of electrostatic potential as a predictor of hydrogen bonds and other weak interactions as opposed to the current geometrical definitions used today.

Electrostatic Potential
The molecular electrostatic potential (MEP) is an integral part in linking the results found in EDD studies such as this to real physical outcomes such as atom reactivities and how electrons move when weak interactions are formed. The MEP maps shown in Figures 11-13 are a visualisation of these properties, and a comparison of these maps for the molecules in their individual and co-crystalline states provides insights into which groups are involved in the formation of weak interactions and the impetus behind them to form. The MEP maps have been mapped onto an isosurface of ρ, which have been plotted on the same scale for comparability. Visual inspection of the maps for (1), (2) and (3) did not yield any unusual features. The most electronegative regions are located near the hetero atoms, and the most electropositive regions were found near the hydrogen atoms. This is expected as most weak interactions found were traditional hydrogen bonds, and these areas correspond to regions where much of the intermolecular bonding occurs in the form of hydrogen bonds. Further analysis of the map of (3) shows excellent complementarity between bond donors and acceptors with them being more electropositive and electronegative, respectively. It was also found that larger differences in electrostatic potential between bond donor/acceptor pairs corresponded to stronger hydrogen bonds, as discussed above. Additionally, the atoms on the periphery of the asymmetric unit with notable electrostatic potential charges were the same atoms involved in intermolecular bonding as discussed previously. The findings from the analysis of the MEP maps in Figures 11-13 validate the future use of electrostatic potential as a predictor of hydrogen bonds and other weak interactions as opposed to the current geometrical definitions used today.

Lattice Energies
The lattice energies of (1), (2) and (3) were determined for the EXP model to assess the effect of the changes in charge distribution during the co-crystallisation process as described above on the overall stability of the system. This has connections to real-world physical outcomes such as the

Lattice Energies
The lattice energies of (1), (2) and (3) were determined for the EXP model to assess the effect of the changes in charge distribution during the co-crystallisation process as described above on the overall stability of the system. This has connections to real-world physical outcomes such as the

Lattice Energies
The lattice energies of (1), (2) and (3) were determined for the EXP model to assess the effect of the changes in charge distribution during the co-crystallisation process as described above on the overall stability of the system. This has connections to real-world physical outcomes such as the Figure 13. Molecular electrostatic potential maps of (4) mapped on an isosurface of ρ.

Lattice Energies
The lattice energies of (1), (2) and (3) were determined for the EXP model to assess the effect of the changes in charge distribution during the co-crystallisation process as described above on the overall stability of the system. This has connections to real-world physical outcomes such as the solubility of the co-crystal complex in solution compared to its individual counterparts. Two methods were used to calculate the lattice energy; the LATEN option in XD2006 [19] based on the method for calculating total interaction energies suggested by Volkov et al. [42,43] and the PIXEL method developed by Gavezzotti [44], which partitions the electron density into 'pixels', which are then used for the calculation. The PIXEL method utilises a theoretical approach where the electron density of the asymmetric unit is calculated and its neighbours generated by symmetry. Thus, changes resulting from intermolecular bonding such as charge redistribution may not be reflected in its results. Refer to Table 4 for the lattice energies calculated for all structures using both methods. The PIXEL calculation for (3) was split into three different pairs: two pairs of the PCM-44BP complex and the two PCM molecules in the middle. The value reported in Table 4 is the sum of the three separate calculations. The LATEN method consistently reports values considerably higher than that of the PIXEL method. This phenomenon has previously been attributed by Spackman [45] to certain intermolecular BCPs not being found even though the Morse equation has been satisfied. This is thought to be due to the extremely flat nature of these critical points. Our results are in accordance with these conclusions with the lattice energies of (1) differing by the smallest amount (10 kJ·mol −1 ) as there is only one molecule in the asymmetric unit, while (2) and (3) both differ by at least 40 kJ·mol −1 and have two and four molecules in the asymmetric unit, respectively.  (2) and (3)  In both methods, the co-crystal was found to be more stable than PCM by approximately 90 kJ·mol −1 . Considering enthalpic and entropic factors associated with solvation, it would be plausible to conclude that a larger net lattice energy also indicates a greater net solubility. The molecular dipole moment (MDM) is an indicator of the spread of charge across a molecule and may also provide insights into the likelihood of solvation of the system. The co-crystal was found to have the largest MDM (12.34D), across the asymmetric unit, while the value for PCM was slightly smaller (10.29D). The large value of the MDM for PCM can be attributed to the strong intermolecular hydrogen bonds it forms, which causes charge to accumulate at the ends of the molecule and depletion in the middle. The value for the SP model (4.52D) lends credence to this conclusion and is in accordance with findings reported by Coppens et al. [46]. The MDM value for (3) is only marginally larger than that of (1). This difference can be attributed to a number of factors including more numerous hydrogen bonds, the nature of the 44BP molecule and the arrangement of the molecules inside the asymmetric unit of (3). The 44BP molecule inherently has a low MDM (3.18D for the asymmetric unit and 2.20 and 2.01D for the individual molecules) due to the inherent symmetry and as a result does not contribute significantly to the MDM. The lack of change in polarisation between the 44BP molecules in (2) and (3) also means it does not contribute significantly to the MDM of (3). Further analysis of the arrangement of molecules in the asymmetric unit of (3) shows that for both PCM and 44BP, the molecules are perpendicular to each other, resulting in a negation of the charge separation due to intermolecular bonding. Evaluation of the MDM of the two PCM-44BP pairs in (3) yielded MDM values of 9.49 and 5.75D, highlighting that each pair individually has some charge separation; however, they seem to negate each other. All the above factors contribute to (3) only having a marginally larger MDM than (1). Evaluation of the MDM for the SP models of (1) and (3) yielded different results with values of 4.52 and 17.22D, respectively. The gas phase nature of the SP calculations does not account for electron density redistribution from intermolecular hydrogen bonding and explains the significantly different MDM values.

Conclusions
Here, we present an EDD study of PCM (1), 44BP (2) and the PCM-44BP (3) co-crystal with the aim of exploring a potential contributing factor towards the poor tableting properties exhibited by Form I PCM, which is currently overcome by co-formulation with gelatine or starch, while co-crystallised complexes perform relatively better. PCM has not been hindered by this limitation as seen through its commercial success; however, it is hoped that the insights gained from studying the weak interactions in the PCM and co-crystal system can be applied to other drugs to improve their physical properties. High resolution X-ray crystallography was used to examine these systems at the electronic level. Topological analysis revealed a total of 4, 8 and 13 hydrogen bonds for (1), (2) and (3), respectively, with π-π interactions also present. Many of the weak interactions found in (1) were also found in (3) and, in both cases, were responsible for the zig-zag packing motif in their structures. 44BP was found to intercalate between the PCM molecules in (3), resulting in a shallower zig-zag packing shape potentially translating into (3), exhibiting better tableting properties, allowing the sheets of PCM and 44BP in the co-crystal to exhibit lateral shear similar to Form II PCM. The weak interactions present in all structures were also reflected visually in the MEP diagrams and the MDM. A larger number of weak interactions was found in (3) compared to (1) and (2), and these interactions ranged in strength from 4.08-84.33 kJ·mol −1 . These interactions allowed the PCM-44BP planes in (3) to be stable; however, the weak π-stack interaction between 44BP molecules allows lateral shear to occur, which potentially gives rise to desirable tabletting properties as described for Form II PCM. When comparing the integrated atomic basin charges of partitions of the PCM molecule, we found that the more numerous weak interactions in (3) resulted in a charge redistribution with the hydroxyl and amide group becoming more negative and the aromatic region of PCM more positive. This resulted in a larger MDM in (3) (µ = 12.34D), indicating an increased polar surface area, and along with the largest lattice energy, provided a good basis to predict how the co-crystal may exhibit more desirable solubility and stability profiles compared to (1).