A digital workflow from crystallographic structure to single crystal particle attributes for predicting the formulation properties of terbutaline sulfate†

A detailed inter-molecular (synthonic) analysis of terbutaline sulfate, an ionic addition salt for inhalation drug formulation, is related to its crystal morphology, and through this to the surface chemistry of the habit faces and surface energy of the whole crystal. Coulombic interactions between the terbutaline cations and sulfate anions contribute 85% of the lattice energy, hydrogen bonding and dispersion interactions contribute 15%. Morphological prediction identifies a plate-like morphology composed of the forms {010}, {100}, {001} and {11̄0} in good agreement with crystals grown from solution. Synthonic modelling of the intermolecular interactions, on a crystal face (hkl)-specific basis, reveals that the surface interactions on the forms with less, relative surface area: {100}, {001} and {11̄0} manifest a greater surface energy, associated, notably, with a greater proportion of polar interactions at these surfaces compared to the morphological dominant {010} surfaces. The predicted total surface energies and their dispersive contributions are found to be in good agreement with those measured at high surface coverage using inverse gas chromatography (IGC). The modelling approach is complementary to IGC, as it enables surface sites of lower energy to be probed, that would require experimentally unachievable surface coverages. The utility of synthonic modelling, in understanding the surface properties of pharmaceutical materials, is highlighted through a workflow-based pathway for digital drug product design encompassing molecule structure, intermolecular packing, crystal morphology, surface energy and formulation properties.


Introduction
Tert-Butyl [2-(3,5-dihydroxyphenyl)-2-hydroxyethyl] ammonium hemisulfate or terbutaline sulfate (TBS) is a β-adrenergic agonist: protonated terbutaline hemisulfate. TBS is present as the active pharmaceutical ingredient (API) in a number of commercial formulations for asthma therapy. 1 The compound has five crystallographic forms 2 including two anhydrous modifications (A and B), one monohydrate, one higher hydrate and one ethanoic acid solvate. 3,4 The stable modification of TBS is an anhydrate B (Fig. 1) and it is this structure which is the subject of the current investigations. Two single-crystal structure determinations are reported for anhydrate B (ZIVKAQ 5 and ZIVKAQ01 (ref. 6)) at 293 K which have similar R-factors (6.23% and 7.30%, respectively). Commercially, TBS is marketed as a racemic compound which decomposes on melting with the melting point being found to lie in the range of 258-260°C. 7 The commercial drug formulation typically contains respirable TBS drug particles that have a particle size of ∼1-5 μm, formulated with an α-lactose monohydrate carrier crystals having a particle size of 70-140 μm. [8][9][10] Due to their small particle size, respirable API particles display highly cohesive and adhesive behaviour readily tending either to form agglomerates 11 or to adhere to the surfaces of excipient materials, manufacturing process equipment or inhaler device components. [12][13][14][15] Agglomerated or adhered drug particles must be re-dispersed effectively upon inhalation in order to be suitable for lung deposition. Hence, a detailed understanding of their surface chemistry and surface energy is important, in terms of determining the magnitude of interparticulate interactions, for designing better formulations e.g. improving content uniformity, aerosolization and inhaled product performance.
Indeed, there have been relatively few studies of the effect of variations in the particle production processes upon surface crystallinity, particle adhesion and cohesion in terms of the inhalation performance of TBS. 10,16,17 The work carried out so far has not provided any fundamental insight as regards the surface properties of TBS from a molecular-scale standpoint, particularly in terms of characterizing the external morphology, surface chemistry and particle surface energy. TBS is an interesting compound of generic relevance in the field of inhalation pharmaceutics. Detailed knowledge, provided by a fundamental intermolecular (synthonic) assessment, can deliver useful baseline data for improved understanding of inhalation-particle surface properties and associated interparticle interactions. 18 Accurate characterization of particle surface-energies can be challenging. [19][20][21] For example, sessile-drop contact angle measurements on the individual facets of macroscopic single-crystals and/or powder compacts via the capillary rise method have not always been found to provide reliable results. 22,23 Surface energy measurements from inverse gas chromatography (IGC) using both polar and apolar probes 17,[24][25][26][27] have, however, been found to be particularly useful in studies of the crystallinity, surface energy and surface properties of particles. Elsewhere, work using IGC has revealed a number of problematic issues such as sample preparation-dependency, large sample sizes (e.g. for micronised materials ∼50-100 mg and for coarse particles more than 2 grams of sample are required typically), 28 timeutilisation and hence high cost. Variability in IGC measurements has been reported previously [29][30][31] and the technique, even with finite dilution analysis, has not always been found to reliably assess the totality of a powder's surface area. The later reflects the technique's tendency to preferentially characterize the higher energy binding sites. 32,33 Hence, if results are to be useful, there is a clear need for a better molecular-scale understanding of both the sources of variability and the impact of surface properties on IGC measurements.
Previous work by Grimsey et al. 34 has shown the utility of molecular modelling in terms of understanding the surface chemistry of drug particles and also, its capability for identifying the interaction sites of molecular probes involved in IGC analysis. In this previous work, the micronized material was examined by making the assumption of particle breakage parallel to planes having the lowest surface attachment energies, i.e. the morphologically dominant crystal surfaces. 35,36 This is despite the possibility that fracture could take place on non-crystal habit surfaces which might in turn have high surface energies and rugosity. 37 Molecular dynamics simulations have also been employed to calculate the surface adsorption energies of the probe molecules, 38 giving a reasonably close match to those determined from IGC. However, these studies were only carried out for the largest crystal surfaces, i.e. those which generally have the lowest surface energies. So far, no studies have adopted the holistic approach of integrating contributions from all the habit faces that comprise the external morphology of the crystal.
Synthonic modeling has been exploited extensively to investigate lattice energy, attachment energy and for the prediction of crystal morphologies using empirical forcefield methods. [39][40][41] The interior and surface properties of crystals can be predicted 42,43 through calculations of the strength and direction of the inter-molecular interactions associated with the bulk (intrinsic synthons) and surface terminated (extrinsic synthons) structures. 41,[44][45][46][47] These approaches have been used to investigate a pharmaceutical salt 46 and also to examine the cohesivity of APIs and excipients. 18 However, acid addition salts (of which TBS is an example) have not, as of yet, been characterised despite their importance as APIs for inhalation drug formulation. Synthonic modelling tools can also predict the surface energy of a crystal and hence provide a clear link to surface energies measured by IGC. 48,49 This modelling approach is also useful in being able to partition the calculated surface energy between the different contributing crystal surface planes (hkl) present on a fully facetted drug particle. Such attributes have obvious importance in terms of aiding the understanding of blending, flow and granulation for both APIs and excipients.
In summary, this paper reports an analysis of the key solid-state synthons, together with the predicted morphology, surface chemistry and surface energy of the morphologically important crystal surfaces of TBS. This analysis is complemented with calculations of solvent probe interactions with the TBS crystal habit surfaces to relate these both to the measured surface energies from IGC and to the surface chemistry that characterizes the crystal morphology.

Materials
The molecular formula of TBS within the asymmetric unit is 2ĳC 12 H 20 NO 3 ] + SO 4 2− with molecular mass, M w = 548.65 g mol −1 . 5 The anhydrate form B of TBS crystallises in a triclinic structure with space group P1 with unit cell parameters a = 9.968 Å, b = 11.207 Å, c = 13.394 Å, α = 100.86°, β = 104.42°, γ = 101.63°with a tri-ionic (2 cation and 1 anion) asymmetric unit 5 (see Fig. 1). In this work, the atomic coordinates were taken from the crystal structure 5 with the lower R factor (refcode: ZIVKAQ). For the experimental work; terbutaline sulfate was purchased from Sigma Aldrich (purity ≥98%) as a white to grey-white crystalline powder for crystallisation experiments. For surface energy measurements, terbutaline sulfate (batch AAEX) was kindly provided by AstraZeneca (Mölndal, Sweden), and was used as supplied.

Crystal growth
As received TBS was purified by recrystallization. TBS single crystals were prepared by cooling crystallisation as follows: 0.32 g of TBS were added to a mixture of 70% water and 30% ethanol (% by mass) to make up 1 g of solution. This was then heated to 70°C and stirred for 30 minutes to ensure complete solute dissolution and then cooled to 5°C and left for 2-3 days until crystallization occurred. An Olympus BX51 microscope, with a UMPlanFl 10×/0.30 objective was used to characterize the size and shape of the resultant crystals.

Surface energy measurements
Specific surface area (SSA BET ) and surface energy (SE) analysis was conducted using a surface energy analyser (iGC-SEA, Surface Measurement Systems Ltd, UK; https://www. surfacemeasurementsystems.com/products/igc-sea/). 0.247 g of terbutaline sulfate was packed into a silanised glass iGC column (internal diameter 4 mm). Prior to any measurements, the loaded column was conditioned using helium carrier gas at 10 scc min −1 for 2 h at 30°C and 0% RH. Methane gas was injected at the start and the end of the experiments for the dead volume calculation. SSA BET was calculated via Brunauer-Emmett-Teller (BET) theory, based on the n-octane adsorption isotherm data and using the peak max parameter. 18 For the surface energy, the columns were equilibrated as mentioned above. Non-polar probes (n-decane, n-nonane, n-octane and n-heptane) and polar probes (chloroform, toluene, ethyl acetate and acetone) were injected into the column at concentrations consistent with a range of surface coverages (n/n m , between 0.5 up to 13% for all probes with the exception of n-decane that reached maximum surface coverage at 4.5%). Based upon an analysis method previously reported, 50,51 the dispersive (γ d , non-polar) and acid-base (γ ab , polar) components of the surface energy were calculated using the Dorris-Gray method and the peak centre of mass parameter. This is based on a linear fit on the 4 alkane probes for surface coverages 0.5-4.5% and based on 3 alkanes for surface coverages above 4.5%. Specific free energies of adsorption (ΔG SP ) of the monopolar probes were obtained using the polarisation approach 52 and the peak centre of mass parameter. All measurements were made in triplicate.

Computational modelling methodology
The computational analysis was carried out utilizing the Cambridge Crystallographic Data Centre (CCDC) materials Mercury, 53 BIOVIA Materials Studio, 54 HABIT98 (ref. 39) and SystSearch 55,56 software. Partial atomic-charges were calculated using the semi-empirical quantum mechanics program MOPAC 57 utilizing the Austin model 1(AM1) approach. The energetic strength and chemical identity of the intermolecular synthons, the lattice energy (E latt or E cr ), and, hence, the slice (E sl ) and attachment (E att ) energies were calculated with HABIT98. 39 This was aided by a detail analysis of the intermolecular interactions (synthons), both within the bulk structure and on the crystal habit faces, which was complemented by molecular-scale visualization of the bulk and surface chemistry. The digital workflow summarized in Fig. 2 provides an overview of the main features of the analysis used.
2.4.1 Calculation of bulk (intrinsic) synthons and lattice energy determination. The lattice energy was calculated using potential parameters from Tripos 5.2, 58 which employs an empirical interatomic potential, Lennard Jones 6-12 potential, see eqn (1), with a scaling factor for hydrogenbonds, together with the Evjen summation method 59 to ensure convergence of the coulombic energy summation.
where A and B are atom-atom parameters, i is an atom in the central molecule, j is an atom in the kth surrounding molecule, N is number of surrounding molecules, n is total number of atoms in the central molecule, n′, is total number of atoms in each of the surrounding molecule, q i and q j are partial, atomic point charges on atoms i and j, and r ij is inter-atomic distance between atoms i and j. The lattice energy was calculated using spherical radii, for cut-off of the energy summation, between 5 and 50 Å, in 5 Å steps, to test and ensure convergence of the lattice energy. The calculation was repeated for each molecule in the asymmetric unit and the results were averaged over all sites in the asymmetric unit and then averaged over both the asymmetric units in the unit cell.
2.4.2 Calculation of surface (extrinsic) synthons and morphological prediction. The pairwise intermolecular interaction energies between a central asymmetric unit in the crystal lattice and all the surrounding asymmetric units, contributing to the lattice energy summation, were partitioned into slice and attachment energy contributions, according to eqn (2). If the centre of coordinates of the remote asymmetric unit was within the region bounded by adjacent Miller planes (hkl), separated by the inter-planar spacing d hkl , then the associated pairwise interaction-energy is regarded as contributing to the slice energy (E hkl sl ) for that crystal surface plane (hkl) orientation. Conversely, if the remote asymmetric unit is centred outside this region the associated interaction energy is taken to contribute to the attachment energy E hkl att .
The attachment energies were normalised to that calculated for the slowest growing face, i.e. with the largest surface area to determine relative growth rates per surface. The Mercury software was used to construct a Wulff plot 60 and, through this, a prediction of the crystal habit based on the attachment energy model was generated. The latter assumes that the relative growth rate (R hkl ) of a given crystal surface (hkl) can be taken, at low supersaturation, to be proportional to E hkl att , 44,61,62 see eqn (3).
The morphological predictions were crossed-validated through experimental studies of the observed crystals as obtained using cooling crystallisation techniques (see section 2.2). The anisotropy factors for the crystal habit surfaces were calculated using eqn (4), where the anisotropy factor, ξ hkl , provides a useful indication regarding the degree of synthon saturation due to surface termination effects. [63][64][65] The anisotropy factor provides a measure of the fraction of intermolecular interactions that are saturated, for a molecule exposed on a particular crystal surface (hkl), when compared to same molecule within the bulk crystal structure.

2.4.3
Calculation of surface energy. The surface energy (SE) of crystal surfaces was calculated from the attachment energy of individual crystal surfaces 62 using eqn (5).
where Z is the number of molecules per unit cell, d hkl is the thickness of the growth step layer, N A is Avogadro's constant, and V is the unit cell volume. The dispersive and polar components of the surface energy were calculated from the surface energy contribution from the van der Waals and coulombic components, respectively. The overall particle surface energy (γ particle ) 49 was calculated using eqn (6), as a surface area weighted average based on the attachment energy.
Here n is the number of forms ({hkl}, i.e. family of faces), γ hkl i is surface energy of individual crystal surfaces, A hkl i is the fractional surface area of the habit face and M hkl i is the multiplicity of the habit face.
2.4.4 Calculation of interaction energies for heterogeneous probe molecules on the crystal surfaces. The binding of heterogeneous probe molecules on the crystal habit surfaces was predicted using systematic grid search methods, 18,66-69 using the SystSearch 18 program. The model is based on a calculation of the rigid body intermolecular interaction energies between selected probe molecules and crystal surfaces. This involves the generation of a three-dimensional grid adjacent to the model of the crystal surface with the interaction energies being calculated by locating the probe molecule at every grid point on and above the surface where the probe is rotated through a grid of points in a rotational Eulerian space. The calculations of the interaction energies were carried out without surface relaxation using the Dreiding II atomic force field, 70 and with partial charges calculated using the MOPAC method. 57 The forcefield contains separate functions for calculating the van der Waals (dispersive), hydrogen bonding (H-bonding) and coulombic interactions and hence the respective contributions from each type of these interactions to the total intermolecular interaction was quantified. A more comprehensive description of this methodology was provided by Ramachandran, V., et al. 18 From this, the strongest interaction (most negative interaction energy) and the mean interaction energy from the histogrammic distribution 'finger-print' was calculated for each crystal habit surface (hkl).

Results and discussion
Results from the molecular and synthonic modelling analysis are described in sections 3.1 and 3.2, the morphological analysis and surface chemistry are presented in section 3.3 and 3.4, the comparison of the calculated surface energies with those measured using IGC is given in section 3.5 and the cross-correlation of the IGC data with intermolecular binding energies is given in section 3.6.

Analysis of intermolecular packing within the bulk crystal structure
The intermolecular packing and crystal chemistry of TBS is summarised in Fig. 3 and Table 1.
The terbutaline cations have an asymmetric centre at carbon atom C7 and C19 for terbutaline cation 1 (TB I) and terbutaline cation 2 (TB II), respectively. The terbutaline cations consist of phenethylamine (C 6 H 5 -CH 2 -CH 2 -NH 2 -) groups in which the nitrogen atoms N1 and N2 are protonated (see Fig. 1). Each terbutaline consists of three alcohol groups (two -OH groups attached the benzene ring and one -OH group attached to the alkyl chain) and one amine -NH 2 group. The asymmetric unit contains two terbutaline cations and one sulfate anion. The two crystallographically inequivalent terbutaline ions (TB I and TB II) exhibit differences in their conformation reflecting different rotations of the tertiary-butyl group about the connecting nitrogen to carbon covalent bond (NH 2 -C-(CH 3 ) 3 groups).
The basic unit of the crystal structure is an inversionrelated terbutaline ion TB I dimer formed between the amine and alcohol groups attached to the alkyl chain. Each terbutaline cations TB I is bonded to one cation TB I, one cation TB II and two sulfate anions (Fig. 3a). Each cation TB II is bonded with one cation TB I and one cation TB II and 3 divalent sulfate anions (Fig. 3b).
In the crystal structure, each sulfate anion forms six hydrogen bonds with five TB cations (Fig. 3c) (two TB I and three TB II, one of the two, TB I, forms 2 hydrogen bonds with the sulfate group through N-H⋯O-S and O-H⋯O-S) utilizing six of the eight potential hydrogen bond acceptors from sulfate group. Each TB cation donates five (from three -OH groups and one -NH 2 group in each TB cation) and accepts two hydrogen bonds (from two -OH groups).
In summary, there is a strong hydrogen bond network within the TBS crystal structure, with the use of all ten the hydrogen bond donors. However, not all the potential hydrogen bonding acceptor sites are used in the crystal structure, with four acceptors from TB I and TB II (two nitrogen atoms in the -NH 2 group and two oxygens in the -OH group) and two acceptors from sulfate group are being unsaturated.

Lattice energy and identification of key bulk intrinsic synthons
Partial atomic-charges, calculated using the semi-empirical quantum mechanics program MOPAC, are provided in the section S1 in the (ESI †). Summation of the intermolecular interactions with distance revealed convergence at 15 Å to a lattice energy of −138.88 kcal mol −1 of asymmetric units with the coulombic interactions accounting for ∼85% of the lattice energy (Fig. 4). To the authors' knowledge, no sublimation enthalpy data has been published for TBS with which to compare to the calculated lattice energy. The summation of the percentage contribution to the lattice energy from the different molecular groups of the asymmetric unit, revealed that, TB I and TB II groups contributed 18.4% and 22.8% respectively, whilst the SO 4 2− group was found to make the greatest contribution ca. 58.8% reflecting its intimate involvement with the strong coulombic anion-cation interactions. Table 2 lists the top nine strongest, attractive intrinsic synthons found from the analysis of the bulk crystal structure and highlights the component of ions in the attractive synthons (with subscript a). All of the most attractive synthons were found to be between the sulphate anionterbutaline cation pairs, indicating that the strongest interactions within crystal lattice are a result of the coulombic forces. Also, due to the symmetrical inequivalence of the terbutaline cation pairs, the synthon types found are quite similar to each other and simply occur between the different symmetrically inequivalent ions. Fig. 5 reveals that the orientation of specific atomic groups, particularly the -OH and -NH 2 functional groups strongly influences the strength of the synthon. Synthons A a , B a , C a and D a (interactions between sulfate with TB I or TB II) have the closest orientation of the oppositely charged SO 4   the SO 4 2− and the NH 2 + groups are closer in synthon L a , synthons I a /K a and J a were found to have a greater attractive interaction due to these synthons forming more significant van der Waals and H-bonding interactions, respectively. The strongest repulsions within the crystal structure were found to be between the sulfate anions (see ESI, † section S2). Due to the symmetrical equivalence of the anions, the four most-repulsive, independent pairs only are illustrated. Table 2 lists that, as expected, the majority of the repulsion energy is due to the coulombic repulsions. These strongest repulsive synthons were found to be between 6.43-9.48 Å apart (S-S). These repulsions were more than compensated for by the strong attractive forces between the anion-cation pairs.

Predicted morphology and identification of key surfaceterminated extrinsic synthons
The predicted morphology, shown in Fig. 6(a) and detailed in Table 3 column 1-6 based on the attachment energy model, shows a plate-like morphology with the {010} being the dominant form, along with the smaller {100}, {001} and {110} surfaces. The predicted morphology agrees well with the observed experimental crystal morphology of TBS crystals grown in a mixture of 70% water and 30% ethanol at 5°C (see Fig. 6(b)).
The percentage surface saturation (anisotropy factors) are listed in Table 3 column 7. As expected, the {010} surface contains the highest degree of saturation, 96.28%, due to fewer strong synthons being broken at the surface. In Table 2 Strength of top nine most attractive (with subscript a) and 4 most repulsive synthons (with subscript r) along with synthon analysis for the {010}, {100}, and {001} and {1−10} crystal surfaces. Distances quoted are the distances between the two centres of gravity of the ions. Note that SF indicates sulfate groups and TBI and TB II are terbutaline cations 1 and 2, respectively; SL is slice and ATT is attachment

Ions involved
Inter-ion distance (Å) Van Table 2 indicating that the major interactions are between the oppositely charged ions. Similar interaction types that are only different because they are between the different terbutaline cations type I and type II all distances are quoted between S in SO 4 2−  In this study the synthons that are predicted to contribute to the attachment energy, and hence growth, of a given crystal surface have been identified and displayed in Table 2. This reveals that stronger intermolecular interactions (synthons) contributing to the attachment energy E att of the crystal faces {hkl} leads to the higher growth rates of {hkl}. The strongest interaction contributing to the attachment energy at the {010} surface was found to be synthon F a (−41.37 kcal mol −1 ) (see Table 2), which is the coulombic interaction between sulfate and amine group of terbutaline cation 2. The strongest interactions contributing to the attachment energy at the {110} surface were found to be E a , F a , and G a (−48.45, −41.37 and −39.00 kcal mol −1 ). These predictions are consistent with the morphological dominance of this face and with the order of the importance of the faces:

Analysis of the surface chemistry of the habit planes
The morphology of a crystal in different growth environments depends, collectively, on the relative growth rates associated with its constituent crystal-habit faces. Ultimately these rates are linked to the respective interaction energies associated with the strongest synthons on these surfaces. The growth morphology in a vacuum of an individual habit face can be rationalised in terms of the respective molecular arrangements, and hence synthon orientation, on these specific faces. Hence, molecular modelling can be utilised to examine the surface chemistry of the {010}, {100}, {001} and {110} surfaces. For the faces listed in Table 3, the extrinsic synthons identified as contributing to the attachment energies were labelled on molecular models of the surface as shown in the section S3 in the ESI. † The {010} form exhibits one -OH group from each terbutaline cation exposed on the surfaces with all the H-bonding being embedded within crystal surfaces. The surfaces of the forms {100} and {001} exhibit two -OH groups one each from each terbutaline cation exposed on the surfaces; however, the {001} surfaces have non-bonded H-bond acceptor groups across the surface which are available to form hydrogen bonds with probe molecules. The {110} morphological form exposes three -OH groups and the protonated -NH 2 -C 3 H 9 moiety on the crystal surface. This indicates that these faces are more polar than the other faces and explains the ranking of the morphological importance of the exposed crystal surfaces in terms of surface area.

Comparison of predicted cumulative particle surface energies with measured IGC data
The surface energies, as a function of surface plane, together with the fractional surface areas of the associated habit plane are given in Table 3, columns 9 and 4, respectively. The predicted cumulative surface energies and the measured surface energies of non-micronized TBS, as measured by IGC, are reported in Table 4 and Fig. 7. The dispersive component of the surface energy measured by IGC ranges between 49.1  The predicted surface energy showed significant differences between the surface energies of individual crystal surfaces (see Table 3, columns 8 and 9). The surface energies of the {110} and {001} forms are greatest whilst the surface energy of the {010} form is smallest. This agrees with the surface chemistry analysis which indicates that the {110} and {001} forms have more polar and non-polar components exposed on the surfaces contributing to the dispersive and polar surface energy, for example, for the {110} faces the -CH 3 -CH 3 -CH 3 and -CH 2 group (alkyl groups) attached to the -NH 2 contribute to dispersive energy components and -NH 3 and three -OH groups contribute to polar surface energy components.
The predicted surface energy, which is based upon an idealised perfect termination of the bulk crystal lattice, gave smaller values than those measured for TBS supplied from AZ (see Table 4) and those measured for micronized and processed TBS reported by Rehman, M. et al. 71 It has been reported that IGC performed in the infinite dilution regime usually yields higher values of dispersive surface energy than other measurement techniques, 72 as well as sometimes providing an unrepresentative measure of the majority of the powder surface area. Infinite dilution IGC, historically performed at 0.03 p/p 0 , has been shown to probe the crystal surface regions with the highest energy, while finite dilution IGC can provide a surface heterogeneity map. 32,[73][74][75] The effectiveness of infinite dilution IGC is likely due to adsorbed molecules having no other adsorbate molecules close enough to allow any adsorbate-adsorbate interactions, hence there is a lack of competition for adsorbate binding sites. The highest energetic sites probed by infinite dilution IGC may not dominate the inter-particulate contact areas, since their relative contribution to the total particle surface area may be exceedingly low. For example, for TBS the {001} faces contribute the highest surface energy (Table 3), but this face only provides 8.69% of the crystal surface area (Table 3).
It was previously reported 33,76 that size reduction, such as micronization, leads to an increase in the measured dispersive surface energy of the drug substance with increasing extent of micronization. The increase in surface energy is probably due to the generation of new, higher energy sites of crystalline disorder due to impaction events or particle fracture. 33 Infinite dilution IGC measurement of surface energy would therefore be expected to be dominated by the latter, unrepresentative high energy sites.
It is significant that computational prediction of surface energy for the highest energy sites (the {001} surfaces: 109 mJ m −2 ) correlated well with the experimental measurements at low surface coverage (103.3 mJ m −2 ), where the highest energy surface sites are probed (both from the standpoint of crystallography and process-induced disorder). It was evident that the surface energy values determined experimentally for the highest surface coverage (which should probe the relatively lower surface energy sites) showed an improved  correlation with the average calculated surface energy when this is weighted with respect to the relative contribution of each crystal face to the total surface area of the crystal (see Table 4).
3.6 Comparison of the predicted adsorption energy for the different polar probes with those measured using IGC The results for the grid-search modelling are given in Fig. 8. These provide histograms of calculated interaction energies for a series of probe vapours (chloroform, toluene, acetone, and ethyl acetate) on TBS crystal surfaces {010}, {100}, {001} and {110}. The strongest interaction energies, the mean interaction energy distribution of a probe molecule on TBS crystal surfaces, together with the total interaction energy of probe molecules on the TBS particles are summarized in Table 5. This study elucidates the anisotropic surface chemistry of the major faces of TBS, identifies that polar solventprobe molecules were found to interact more strongly with the side faces {100}, {110} and {001} and least with the large {010} surface. This illustrates that the relative contribution of a crystallographic face to its total surface area plays an important role in the preferential binding site for the probe molecules. Indeed, the modelling studies indicated that the surface roughness at the molecular-level can have a major influence on the probe molecule adsorption energies, which may impact the actual process of adsorption in experimental determinations. This increased surface roughness can be expected to provide a favourable docking site for adsorption in a surface 'valley', for the probe molecules. Hence, it has been shown that the higher interaction energy sites for the {100}, {001} and {110} surfaces tend to be located in the surface 'valleys'. Particularly for the {110} surface, the probe molecules prefer to bind in a channel present at the crystal surfaces, which is also the case for ethyl acetate, chloroform and toluene probe molecules, resulting in high interaction energies at this surface. Further details for the preferred binding position of the probe molecules on the TBS crystal surfaces are presented in the ESI † (see section S4).
The data analysis also elucidated that for toluene and chloroform, the solvent probe interactions with all the surfaces are found to be almost completely from dispersive force contributions and coulombic interactions (for chloroform), without a hydrogen bonding contribution. However, it should be noted that in a real solvent system, chloroform is known to form weak hydrogen bonds, the Dreiding forcefield used in these calculations does not consider chlorine as a hydrogen bond acceptor atom and hence no hydrogen bonding contribution would be expected. Similarly, toluene as a non-polar solvent molecule, is incapable of forming hydrogen bonds so would also be expected to have zero contribution to the overall interaction energy. In contrast, for ethyl acetate and acetone, the interaction energies are found to have had a significant H-bonding contribution and, additionally, the coulombic energy contribution is significant, particularly for acetone.
The measured specific surface energies ΔG SP (kJ mol −1 ) using chloroform, toluene, ethyl acetate, and acetone are given in Fig. 9.
Despite the poorer agreement, a trend of decreasing measured surface energy was apparent as surface coverage was increased, probing the lower surface energy sites (see Fig. 7b and 9b). The trend for surface energy in Fig. 9 indicates that the surface energy distribution of TBS is quite heterogeneous. During the finite dilution IGC experiments, the surface energy distribution was calculated via the retention volumes of the probes at a specific surface coverage in order to allow for the differences between the various probes. However, this is also a limitation of the technique as the maximum injectable amountand as such the maximum surface coverage a selected probe can achieveis determined by the probe's cross-sectional area. In this present work, it was not possible to measure the dispersive component of TBS surface energy accurately above surface coverages of 13% as that would reduce the number of available alkane probes to two (octane and heptane), insufficient for accurate linear regression. Further analysis of the adsorption isotherms of the probes (see Fig. 7a and 9a, and section S5 in the ESI † for the adsorption isotherms) showed that even at the surface coverage of 13%, the partial pressure in the column did not exceed 0.06 P/P 0 Table 5 Strongest (best binding position), mean and total interaction energies of a probe molecule on TBS crystal surface alongside the surface free energies of adsorption as measured from IGC at 0.5% and 13% of surface coverage  which indicates that a full monolayer coverage has not been reached and thus only a fraction of the heterogeneous surface was probed. Consequently, this would suggest that the experimental surface energy measurements by IGC did not achieve sufficient surface coverage to probe the lowest surface energy sites in the same way as have been predicted in the molecular modelling-based grid search approach.

Conclusions
In this work, molecular and synthonic modelling tools have been successfully applied to link dominant synthon contributions, the percentage surface saturation, the surface area, surface chemistry analysis (functional groups exposed on crystal surfaces, extrinsic synthons) surface energy and physical properties. Synthonic modeling revealed the extent of the importance of coulombic interactions between the salt pairs to stabilize the lattice energy, compared with small, free base, molecule crystal structures. Attachment energy morphology predictions produced a plate-like morphology with the {010} being the dominant form, along with the smaller {100}, {110}, {101} and {001} forms, which is consistent with the crystal morphology observed from solution grown crystals. The synthon analysis explains why the relative morphological importance to the crystal morphology of forms is in the order {110} < {001} < {100} < {010}. The faster growing {110} surfaces were predicted to have many more of the strong synthons notably E a , F a and G a contributing to their attachment energy. The calculated dispersive and total surface energy had lower values than the measured data reported in this study and tended to correlate better at a surface coverage of 13% compared to 0.5% surface coverage. This might be due to the experimental surface energy measurements by IGC in this study which did not achieve sufficient surface coverage to probe all the lowest surface-energy sites, combined with the finding that the surface energy distribution of TBS across its different crystal habit faces is quite heterogeneous. From the calculated adsorption energy for the different polar probes on TBS crystals, the mean interaction energy was found to provide a better correlation with the measured surface free energies of adsorption (−ΔG SP ).
The calculated surface energy analysis was found to be very helpful for assessing and interpreting the measured IGC data, notably through its ability to both partition surface energies between the different morphological forms as well as probe lower energy binding sites not routinely accessible using IGC, for designing formulation processes for inhalation drugs. This is due to the importance of surface properties, their link to drug particulate-interactions and the associated, intrinsic, powder properties needed for achieving effective drug aerosolization, efficacy, lung delivery and therapeutic success. This suggests that the type of analysis presented here might have more general applications in the formulation design of other APIs in order to understand, in-detail, the fundamental material behavior at the molecular scale. Energy of intermolecular interactions found within one d-spacing on the (hkl) crystallographic plane Attachment energy Energy of intermolecular interactions formed when a crystal slab one d-spacing thick, defined by (hkl) plane, is inserted into a bulk crystal structure Anisotropy factor

List of symbols
The degree of saturation of a molecule exposed at a cleaved crystal surface (hkl), in comparison to the same molecule fully saturated in the bulk structure

Conflicts of interest
The authors declare no competing finance interest.