Crystallization: Its Mechanisms and Pharmaceutical Applications

Abstract

The crystallization of small-molecule drugs plays an important role in the pharmaceutical industry. Since many downstream industrial processes are heavily influenced by the crystalline properties of a drug, that is, crystal shape, size distribution, and polymorphic form, control over the crystallization process can facilitate manufacturing and testing. However, before the crystallization process can be controlled, an understanding of its underlying mechanisms is required. In this chapter, we will look at the thermodynamic driving force behind crystallization and how crystal nucleation and growth rates can be used to control the properties of the resulting crystals. Throughout the chapter, we give examples of how these control approaches can be applied in pharmaceutical research and industry to obtain crystals with desired properties. We then finish this chapter with a look at crystallization from the amorphous state, which differs from crystallization from solution and is a relevant topic in pharmaceutical sciences, since the preparation of an amorphous solid is a popular approach to enhancing the solubility of a drug.

1. Introduction

Crystallization plays an important role in the manufacture and purification of small-molecule active pharmaceutical ingredients (APIs). It is estimated that between 70 and 80% of all small-molecule, APIs have at least one crystallization step in their manufacturing processes [1, 2]. To facilitate other downstream processes, such as filtration, drying, dissolution testing, and formulation, it is often desirable to be able to consistently produce crystals with specific properties, such as crystal size distribution, crystal shape (habit), and polymorphic form. Such control over the crystallization process requires accurate descriptions of crystal nucleation and growth kinetics, as well as solubility, breakage, and agglomeration data [3]. Crystallization control has gained even more interest since the release of the United States Food and Drug Administration’s (FDA) Process Analytical Technology (PAT) framework, which aims to improve efficiency in pharmaceutical development, manufacturing, and quality assurance through innovative process development, analysis, and control [4].

Designing and implementing controlled crystallization processes can be time-consuming, and usually involve trial-and-error, model-based, or model-free (direct-design) approaches [5]. Accurate kinetic descriptions of crystal nucleation and growth form the backbone of the model-based approach which, although more complex, has the potential to alleviate the labor-intensiveness and facilitate the optimization of the traditional trial-and-error approach [6]. The models used in the model-based approach can even be used together, to more accurately describe the underlying crystallization mechanism(s). For example, in a study by Quilló and coworkers [6], three different crystal growth models were fit to experimental data, to find the most likely crystallization mechanism. They found that the underlying crystallization mechanism could best be described by an additive combination of simultaneous Birth-and-Spread (B+S) and Burton-Cabrera-Frank (BCF) models, which will be discussed later in this chapter. An example of the model-free approach is direct nucleation control (DNC), which attempts to directly control the number of nuclei present in a system. Using DNC, Abu Bakar and coworkers [5] were able to produce glycine crystals from water–ethanol mixtures that were larger than those obtained from uncontrolled crystallization.

To implement the above-mentioned approaches, the investigator needs to be able to monitor the crystallization process. The instrumental technique used will depend on what aspect of the process one needs to monitor. For the model-based approach, one would typically be interested in monitoring the supersaturation decay during crystallization. This can be done using techniques, such as attenuated total reflectance-Fourier transform infrared (ATR-FTIR) or near infrared (NIR) [7]. If the goal is to monitor the evolution of crystal size distribution over time, for example in the DNC study mentioned above, a technique like focused beam reflectance measurement (FBRM) can be used [5, 7]. In most investigations, the above-mentioned instrumental techniques are used concurrently, to gain as much information about the crystallization process as possible.

The crystallization processes and monitoring techniques mentioned above are concerned with crystallization from solution, where a reduction in temperature and/or the addition of an antisolvent is used to create supersaturated conditions. Under supersaturated conditions, the solute molecules can come together to form small masses, or nuclei, through a process called nucleation. The addition of more solute molecules to these nuclei is called crystal growth. There is also another type of crystallization that is of particular importance in pharmaceutical sciences, and that is crystallization from the amorphous state.

In this chapter, we will discuss the theoretical background of crystal nucleation and growth, look at examples of how these theoretical models can be implemented practically and give examples of applications in pharmaceutical research. We will also look at crystallization from the amorphous state and discuss techniques that can be used to delay it.

2. The driving force behind crystallization

The thermodynamic driving force behind crystal nucleation and growth rates is a chemical potential difference (Eq. (1)):

$${\mu }_i-{\mu }_i^{\mathit{sat}}=\mathit{RT}ln\left(\frac{a_i}{a_i^{\mathit{sat}}}\right)$$E1

where ${\mu }_i$ and $a_i$ are the chemical potential and thermodynamic activity of solute i, respectively, and the superscript “sat” represents the property at supersaturated conditions, R is the universal gas constant and T is the absolute temperature [6, 8]. Because the activity of a solute is hard to measure experimentally, it is usually substituted with the term ${\gamma }_i{x}_i$, where ${\gamma }_i$ is the activity coefficient and $x_i$ the mole fraction of solute i. Since mole fractions can be measured experimentally, Eq. (1) is often rearranged and expressed in terms of a supersaturation ($\sigma$) driving force [8], presented here in Eq. (2):

$$ln\left(\sigma \right)=ln\left(\frac{{\gamma }_i{x}_i}{{\gamma }_i^{\mathit{sat}}{x}_i^{\mathit{sat}}}\right)=ln\left(\frac{a_i}{a_i^{\mathit{sat}}}\right)=\frac{{\mu }_i-{\mu }_i^{\mathit{sat}}}{\mathit{RT}}$$E2

The value of $x_i^{\mathit{sat}}$ can sometimes be determined experimentally from equilibrium solubility studies. However, if the solubility extrema cannot be determined experimentally, it can be estimated from a variety of drug solubility models [9], for example, the Apelblat model [10], or the van’t Hoff–Jouyban–Acree (VH–JA) model [6].

Because the activity coefficients still need to be estimated, many studies make the simplifying assumption of setting the activity coefficients in Eq. (2) equal to 1, thereby expressing the supersaturation driving force only as a mole fraction $\left(x_i/x_i^{\mathit{sat}}\right)$, or as the concentration ratios $\left(C_i/C_i^{\mathit{sat}}\right)$, or $\left(S-1\right)$ where $S=\left(C_i/C_i^{\mathit{sat}}\right)$ [11, 12, 13]. Although these simplifications of the supersaturation driving force have been used successfully, they do not necessarily generalize to all solutes and solvents, and one might need to find a way to estimate the activity coefficient of the solute. Fortunately, there is a way to estimate the denominator$\left({\gamma }_i^{\mathit{sat}}{x}_i^{\mathit{sat}}\right)$ in Eq. (2) in one go.

When dealing with phase equilibria, the traditional thermodynamic reference is a supercooled melt of the pure solute compound. Using this reference, the activity of a solute can be obtained from the enthalpy of fusion at the temperature of interest. In its most general form, this is done through Eq. (3) [14]:

$$ln\left(a_i^{\mathit{sat}}\right)=ln\left({\gamma }_i^{\mathit{sat}}{x}_i^{\mathit{sat}}\right)=\frac{∆H_f\left(T_m\right)}{R}\left(\frac{1}{T_m}-\frac{1}{T}\right)-\frac{1}{\mathit{RT}}{\int }_{T_m}^T∆C_p\mathit{dT}+\frac{1}{R}{\int }_{T_m}^T\frac{∆C_p}{T}\mathit{dT}$$E3

where $∆H_f\left(T_m\right)$ is the enthalpy of fusion at the melting temperature, $T_m$, $T$ is the temperature of interest and $∆C_p$ is the difference in heat capacity between the supercooled melt, $C_p\left(l\right)$ and the solute’s solid-state $C_p\left(s\right)$, presented here in Eq. (4):

$$∆C_p=C_p\left(l\right)-C_p\left(s\right)$$E4

A practical problem with Eq. (3) is that $∆C_p$ has to be integrated down from the melting temperature to the temperature of interest, but far below the melting temperature the thermodynamic properties of the supercooled melt are not experimentally accessible [14]. Looking at the right-hand side of Eq. (3), $∆C_p$ is the only physical property that is difficult to obtain experimentally. This has led to several simplifying assumptions regarding Eq. (3), all concerned with how to handle $∆C_p$.

A common assumption is to completely ignore $∆C_p$ [14], resulting in Eq. (5):

$$ln\left({\gamma }_i^{\mathit{sat}}{x}_i^{\mathit{sat}}\right)=\frac{∆H_f\left(T_m\right)}{R}\left(\frac{1}{T_m}-\frac{1}{T}\right)$$E5

However, when working with temperature ranges normally used in pharmaceutical processes, some approximations obtained from Eq. (5) have shown to be inaccurate [15, 16].

Another common assumption is that $∆C_p$ is constant and can be approximated by the entropy of fusion at the melting temperature $∆S_f\left(T_m\right)$, yielding Eq. (6) [14, 17]:

$$ln\left({\gamma }_i^{\mathit{sat}}{x}_i^{\mathit{sat}}\right)=\frac{∆H_f\left(T_m\right)}{R{T}_m}ln\left(\frac{T}{T_m}\right)$$E6

If we can experimentally determine the isobaric heat capacity of the solid, $C_p\left(s\right)$, at different temperatures below the melting temperature, as well as the isobaric heat capacity of the melt, $C_p\left(l\right)$, the assumption can be made that $∆C_p$ is constant and equal to its value at the melting temperature, $∆C_p\left(T_m\right)$, as illustrated in Figure 1 [14]. Provided that the compound in question does not decompose upon melting, this assumption gives Eq. (7):

Figure 1.

Illustration of the experimental heat capacities (solid lines) and extrapolated heat capacities (dashed lines) used to approximate the heat capacity term in Eq. (3). Reproduced from ref. [14] with permission from Elsevier.

$$ln\left({\gamma }_i^{\mathit{sat}}{x}_i^{\mathit{sat}}\right)=\frac{∆H_f\left(T_m\right)}{R}\left(\frac{1}{T_m}-\frac{1}{T}\right)-\frac{∆C_p\left(T_m\right)}{R}\left(ln\left(\frac{T_m}{T}\right)-\frac{T_m}{T}+1\right)$$E7

If we have enough isobaric heat capacity data above the melting temperature, we can extrapolate down from $C_p\left(l\right)$ to the temperature of interest, see Figure 1, and use Eq. (3) in its general form [14]. In this case, we can use the linear dependency of $∆C_p$ on temperature to rewrite Eq. (4) as:

$$∆C_p=C_p\left(l\right)-C_p\left(s\right)=q+r\left(T_m-T\right)$$E8

where $q$ and $r$ are regression parameters, obtained from extrapolating down from $C_p\left(l\right)$ to temperatures of interest, calculating the difference in heat capacity between the extrapolated data and $C_p\left(s\right)$, and plotting these differences in heat capacity, $∆C_p\left(T\right)$, against $\left(T_m-T\right)$. Notice from the above that the value of $q$ corresponds to $∆C_p\left(T_m\right)$ in Eq. (7). Once we have estimates of $q$ and $r$, we can plug them into Eq. (3) and solve the integrals to give the most comprehensive estimate of $\left({\gamma }_i^{\mathit{sat}}{x}_i^{\mathit{sat}}\right)$ in the form of Eq. (9):

$$ln\left({\gamma }_i^{\mathit{sat}}{x}_i^{\mathit{sat}}\right)=\begin{array}{l}\frac{∆H_f\left(T_m\right)}{R}\left(\frac{1}{T_m}-\frac{1}{T}\right)-\frac{q}{R}\left(ln\left(\frac{T_m}{T}\right)-\frac{T_m}{T}+1\right)\\ -\frac{r}{R}\left(T_{m}ln\left(\frac{T_m}{T}\right)-\frac{T_m^2}{2T}+\frac{T}{2}\right)\end{array}$$E9

Once the values of $\left({\gamma }_i^{\mathit{sat}}{x}_i^{\mathit{sat}}\right)$ at different temperatures have been properly estimated, the values of ${\gamma }_i$ can be calculated from experimental solubility data.

3. Nucleation

In the introduction, we mentioned that crystallization consists of two processes, namely nucleation and growth. Depending on the rates of these two processes, the molecules that make up the crystal may pack differently, giving rise to different crystal structures, or polymorphic forms, of the same compound.

Nucleation is the first step in the crystallization process, and it consists of two mechanisms, namely homogeneous and heterogeneous nucleation. Homogeneous nucleation is triggered by spontaneous fluctuations in the density of the liquid, while heterogeneous nucleation is triggered by contact of the liquid with a foreign solid surface, like an impurity or metastable polymorph [7]. Since most industrial crystallization processes involve heterogeneous nucleation, we will focus mainly on this mechanism, which can be expressed mathematically as Eq. (10):

$$J=N_{0}v\mathit{exp}\left(\frac{-∆G^{\ast }\mathrm{\Phi }}{\mathit{kT}}\right)$$E10

where$J$ is the number of nuclei formed per unit time per unit volume, $N_0$ is the number of solute molecules per unit volume, $v$ is the frequency of molecular transport at the nucleus-liquid interface, $\mathrm{\Phi }$ is the heterogeneous nucleation factor, a function of the contact angle between the nuclei and a foreign surface in the solution with values ranging from 0 to 1, and $∆G^{\ast }$ is the free energy barrier to nucleation of a sphere [7, 18], defined by:

$$∆G^{\ast }=\frac{16\pi {\upsilon }^2{{\gamma }_{\mathit{sl}}}^3}{3{\left(\mathit{kT}\right)}^2{\left(ln\sigma \right)}^2}$$E11

where $\upsilon$ is the molecular volume of the solute, ${\gamma }_{\mathit{sl}}$ is the solid–liquid interfacial tension per unit area, $k$ is the Boltzmann constant and $T$ and $\sigma$ have the same meanings as before [18]. The equation for the rate of homogenous nucleation is similar to Eq. (10) but does not contain the heterogeneous nucleation factor, $\mathrm{\Phi }$, and has a different pre-exponential factor [18].

Looking at Eq. (10), we see some interesting mechanistic features of the nucleation step. First off, we can expect a higher number of nuclei to form from solvents in which the solute is more soluble because these will give higher values of $N_0$, which is a concentration term. Concretely, since the pre-exponential term, $N_{0}v$ is an estimate of the probability of intermolecular collisions, and the term $v$ is mainly determined by the degree of agitation which can be controlled to be constant between different experiments [18], higher solubility will lead to higher values of $N_0$, which will increase the likelihood of nuclei forming. The exponential term itself is a negative exponential, it will decay and asymptotically approach zero for large values in the exponential. In other words, we can expect fewer nuclei to form per unit volume of solute if the free energy barrier to nucleation, $∆G^{\ast }$, is large. The free energy barrier to nucleation is itself dependent on the interfacial tension, ${\gamma }_{\mathit{sl}}$, expressed as:

$${\gamma }_{\mathit{sl}}=0.414\mathit{kT}{\left(C_s{N}_A\right)}^{2/3}\left(ln\left(\frac{C_s}{C_{\mathit{eq}}}\right)\right)$$E12

where $C_s$ is the ratio of the density of the solute to its molar mass, $N_A$ is Avogadro’s number, $C_{\mathit{eq}}$ is the equilibrium solubility of the solute and $k$ and $T$ have the same meanings as above [19]. From Eq. (12), we see that${\gamma }_{\mathit{sl}}$ is inversely proportional to the logarithm of the equilibrium solubility, suggesting that higher values of $C_{\mathit{eq}}$ will result in lower values of ${\gamma }_{\mathit{sl}}$. Putting it all together, we see that under conditions of higher solubility, $C_{\mathit{eq}}$, the interfacial tension, ${\gamma }_{\mathit{sl}}$, will be lower, lowering the energy barrier to nucleation, $∆G^{\ast }$, and thereby increasing the number of nuclei formed, $J$.

In the introduction, we mentioned that there are model-free approaches to controlling crystallization and that DNC is an example of such an approach. In short, the DNC approach attempts to control the number of nuclei present by adding solvent or antisolvent, increasing or decreasing the temperature, or a combination of both, to manipulate the solute’s solubility. Abu Bakar and coworkers employed DNC and were able to produce progressively larger glycine crystals if they lowered the number of newly formed particles [5].

Another example of where the nucleation step of crystallization was used as an intervention point to control the crystal size distribution, is in the work of Fujiwara and coworkers [20]. Paracetamol (acetaminophen), like many other small-molecule APIs, tends to form agglomerates during crystallization, especially if the crystals are small (< 100 μm). This can lead to formulation-related problems further downstream. Fujiwara and coworkers used a solubility curve and determined the metastable limit, that is, the degree of supersaturation that corresponds to spontaneous nucleation. They found that if the degree of supersaturation in a seeded batch crystallization process can be temperature controlled during cooling to stay above saturation but below the metastable limit, as illustrated in Figure 2, larger paracetamol crystals with negligible nucleation and agglomeration can be obtained [20].

Figure 2.

Illustration of how the degree of supersaturation (Run 2) was controlled to fall between the metastable limit (Meta) and the solubility curve (Csat) during cooling. Reproduced from ref. [20] with permission from the American Chemical Society.

4. Crystal growth mechanisms

Once the nucleation step of crystallization is underway, more solute molecules can become incorporated into these nuclei and crystal growth can begin. Generally, crystal growth is considered to consist of two steps in series, namely volume diffusion and surface integration. Volume diffusion consists of the diffusion of solute molecules from the solution, through the boundary layer surrounding the crystal, to its surface. Surface integration can be split into three subprocesses, starting with the desolvation of a solute molecule unto the surface of the crystal, followed by the transfer of this molecule from its point of arrival across the crystal’s surface. The latter is referred to as surface diffusion. Finally, the desolvated molecule reaches an energetically favorable location (kink site) on the crystal’s surface and becomes incorporated into the crystal lattice [13]. The steps in crystal growth described above are illustrated in Figure 3.

Figure 3.

A schematic representation of crystal growth. The steps are (1) the diffusion of a molecule from the solution to the crystal surface, (2) desolvation of the solute molecule unto the crystal surface, (3) transfer across the surface, and (4) incorporation of the molecule into the lattice. The numbering represents—(a) the terrace, (b) the step, and (c) the kink site. Reproduced from ref. [21] with permission from Springer Nature.

In this section, we will discuss the different crystal growth mechanisms, their theoretical basis, and examples of how they are used in pharmaceutical sciences. The topic of surface energy models is quite extensive and will not be covered in detail in this section. However, it plays an important role in how different crystal growth rates can lead to different crystal habits and is, therefore, a relevant topic in solid-state pharmaceutical investigations, as it influences downstream processes, such as dissolution, powder flow, milling, granulation, and compaction. Therefore, we will start this section with a discussion on the Gibbs–Curie–Wulff theorem and then look at a popular surface energy model, namely the Bravais–Friedel–Donnay–Harker (BFDH) model.

Also, in this section, we will not discuss adsorption layer or diffusion-reaction models but will instead focus on the more popular power law, birth-and-spread, and screw dislocation models, and show how information regarding the rate-limiting steps in crystal growth can be inferred from the power-law model.

4.1 Surface energy models

The Gibbs–Curie–Wulff theorem originated when Gibbs postulated that, at equilibrium, a crystal should take on a form such that the product of its surface area and the surface-free energy is minimized. Later, Curie proposed that there exists a direct proportionality between the normal growth rates of crystal faces and surface-free energy. Wulff later stated that, at equilibrium, there exists a central point (the Wulff point) within a crystal such that the distances of the crystal faces from this point are proportional to the specific surface-free energies of those faces [22].

Li and coworkers [22] investigated whether organic crystals, with their diverse molecular shapes, noncovalent interactions, and hydrogen bonding, would still follow the Gibbs–Curie–Wulff theorem. Using a heteroacene molecule as crystal former, they found that not only did the shape and shape evolution of their equilibrium crystals follow the Gibbs–Curie–Wulff theorem, but that the packing of the molecules inside the crystals also reflected the principle of minimizing surface-free energy.

A crystal lattice can be thought of as a repeating 3D pattern consisting of a set of points such that each point has identical surroundings. A lattice has three spatial dimensions ($a,b,c$), also referred to as the lattice constants. In crystallography, Miller indices are used to express a crystallographic plane in terms of a 3D lattice using the notation ($\mathit{hkl}$). The physical meaning of Miller indices ($\mathit{hkl}$) is a set of parallel crystal facets with equivalent interplanar spacing [23]. The interplanar spacing is defined as the distance between two adjacent parallel planes (facets) with the same Miller indices, denoted as $d_{\mathit{hkl}}$. So, for a cubic lattice system, where the lattice constant is simply $a$, the interplanar spacing is defined as:

$$d_{\mathit{hkl}}=\frac{a}{\sqrt{h^2+k^2+l^2}}$$E13

And for a more complex tetragonal system [23, 24]:

$$d_{\mathit{hkl}}=\frac{1}{\sqrt{{\left(\frac{h}{a}\right)}^2+{\left(\frac{k}{b}\right)}^2+{\left(\frac{l}{c}\right)}^2}}$$E14

The BFDH model attempts to predict the crystal habit based on the growth rates of the crystal facets. Concretely, the BFDH model states that the growth rate of a facet is inversely proportional to that facet’s interplanar spacing [21, 25]. Unlike the other crystal growth models that will be discussed later in this chapter, the BFDH model does not have a mathematical representation. Instead, predictions of crystal morphology based on the BFDH model are obtained from computer modeling software, such as Mercury 2020.1, Cerius², or Materials Studio.

Since the BFDH model does not explicitly take into account factors, such as molecular interactions but is instead concerned with the geometric aspects of the crystal, there have been instances where predictions made from the BFDH model were not in agreement with experimental observations. For example, while using the BFDH model to predict the morphologies of two carvedilol polymorphs, Prado and coworkers found that their BFDH model predicted thicker crystals of form II than were experimentally obtained and the model also predicted more faces for form III [26]. Similarly, Nichols and Frampton found that the BFDH model predicted morphologies for paracetamol polymorphic forms I and II that differed significantly from those observed experimentally [27]. Instead of using the BFDH model to predict morphologies, Turner and coworkers [28] used it to identify the important morphological faces of lovastatin crystals obtained from different solvents. They found that the calculated attachment energies of these important faces correlated reasonably well with the observed crystal morphologies.

4.2 Power-law model

The first mathematical model of crystal growth we will discuss is the power-law, or empirical, model. As the name suggests it is an empirical model, raised to some power, and as such it does not offer much insight into the actual growth mechanism. However, it is widely regarded for its flexibility and for being less complex and easier to use than the birth-and-spread or screw dislocation models.

Concretely, the power-law model expresses the growth rate $\left(G_{\mathit{PWR}}\right)$ in terms of the degree of supersaturation as:

$$G_{\mathit{PWR}}=k_g{\left(ln\left(\sigma \right)\right)}^g$$E15

where the kinetics constant, $k_g$, is expected to follow Arrhenius temperature dependence, and is defined as:

$$k_g=k_{g0}\mathit{exp}\left(\frac{-E_a}{\mathit{RT}}\right)$$E16

and where $g$ is the temperature-independent growth order parameter, $k_{g0}$ is a pre-exponential factor, $E_a$ is the activation energy and $\sigma$, $R$, and $T$ have the same meanings as before [6].

We mentioned that the power-law model does not explicitly offer insights into the specific crystal growth mechanism. However, recently certain consistencies between the values of $g$ and $E_a$ and growth mechanisms have come to the fore. For instance, there is now evidence that for values of $E_a$ ranging from 10–20 kJ/mol the rate-limiting step for crystal growth is volume diffusion, and for values ranging from 40–100 kJ/mol, surface integration [13, 17]. For values of $g$ higher than unity, the rate-limiting step is also considered to be surface integration [13].

We also mentioned the flexibility of the power-law model. In a paper concerning the crystal growth kinetics of piracetam polymorphs, Soto and Rasmuson [13] substituted the supersaturation driving force in Eq. (15) with a different driving force, namely a mass balance expression based on Haüy’s law, which states that one can assume a crystal’s shape to be constant if nucleation, agglomeration, breakage and growth rate dispersion can all be assumed to be negligible. The driving force can then be expressed as Eq. (17):

$$ln\sigma =\left(∆C_0-\left[{\left(\frac{\overline{L}}{{\overline{L}}_0}\right)}^3-1\right]\frac{W_0}{M}\right)/C_i^{\mathit{sat}}$$E17

where $∆C_0=C_i-C_i^{\mathit{sat}}$, $\overline{L}$ is the crystal length at any instant, ${\overline{L}}_0$ is the mean initial size of the seed crystals, $W_0$ is the mass of the seed crystals, $M$ is the mass of solvent and $C_i$ and $C_i^{\mathit{sat}}$ have the same meanings as before [13]. Substituting back into Eq. (15) and combining with Eq. (16), they modeled the crystal growth rates using Eq. (18):

$$G_{\mathit{PWR}}=\frac{\mathit{dL}}{\mathit{dt}}=k_{g0}\mathit{exp}\left(\frac{-E_a}{\mathit{RT}}\right){\left[\left(∆C_0-\left[{\left(\frac{\overline{L}}{{\overline{L}}_0}\right)}^3-1\right]\frac{W_0}{M}\right)/C_i^{\mathit{sat}}\right]}^g$$E18

where $\mathit{dL}/\mathit{dt}$ is the rate of change of a crystal’s characteristic linear dimension with time and all the other terms have the same meanings as before [13]. Using Eq. (18), they were able to accurately model the growth rates of their piracetam polymorphs and from the values of $g$ and $E_a$, they were able to determine that surface integration was the rate-limiting step in their crystals’ growth.

4.3 Birth-and-spread model

The birth-and-spread model (B+S) models surface, or two-dimensional, nucleation which is believed to occur at intermediate levels of supersaturation. Under such conditions, surface nuclei can form at the edges, corners, and even faces of crystals, creating new growth steps [13]. These surface nuclei then spread to create a new crystalline layer. The B+S model is presented here in Eq. (19)–(22):

$$G_{B+S}={\tau }_1{\left(\sigma -1\right)}^{2/3}{\left(ln\left(\sigma \right)\right)}^{1/6}\mathit{exp}\left(\frac{-{\tau }_2}{T^{2}ln\left(\sigma \right)}\right)$$E19

$${\tau }_1={\left(\frac{16}{\pi }\right)}^{1/3}{h}^{1/6}{D}_{\mathit{\text{surf}}}{\left({\beta }^{\prime }\frac{{\mathrm{\Gamma }}^{\ast }}{x_s}\right)}^{2/3}{\left(V_m\mathrm{\Gamma }{N}_A\right)}^{5/6}$$E20

$${\tau }_2=\frac{\pi }{3}{V}_{m}h{\left(\frac{{\gamma }_{\mathit{sl}}}{k}\right)}^2$$E21

$$D_{\mathit{\text{surf}}}=A_{\mathit{\text{surf}}}\mathit{exp}\left(\frac{-E_{a,\mathit{\text{surf}}}}{\mathit{RT}}\right)$$E22

where ${\tau }_1$ and ${\tau }_2$ are lumped model parameters, $h$ is the step height, ${\beta }^{\prime }$ is a correction factor $\left({\beta }^{\prime }\le 1\right)$ , $\mathrm{\Gamma }$ is the molecular adsorption coverage of the solute and ${\mathrm{\Gamma }}^{\ast }$ is that property at equilibrium, $x_s$ is the mean displacement of the adsorbed units over the surface, $V_m$ is the molecular volume of the solute, $D_{\mathit{\text{surf}}}$ is the surface diffusion coefficient, $A_{\mathit{\text{surf}}}$ is a pre-exponential factor, $E_{a,\mathit{\text{surf}}}$ is the activation energy for surface diffusion and the rest of the parameters have the same meanings as before [6].

Despite being a very comprehensive model, there is a specific inconsistency between predicted data and experimental observation that should be pointed out, and that is that once growth from a surface nucleation site has spread all the way to the edge of the crystal face, forming a new crystalline layer, any further growth would require adsorption to a smooth surface. Since this kind of adsorption is energetically less favored than binding to dislocations, it would require high levels of supersaturation to overcome the energy barrier. However, experimental data have shown that crystal growth can occur at much lower degrees of supersaturation [8]. To address this inconsistency, a screw dislocation model can be used, which assumes that screw dislocations on a crystal surface serve as self-perpetuating growth sites.

4.4 Screw dislocation model

The Burton–Cabrera–Frank (BCF) model [29] is a screw dislocation model, that assumes that crystal growth stems from screw dislocations on a crystal surface and that growth from these ledges leads to more ledges in a self-perpetuating fashion, thereby attempting to explain how crystal growth can be obtained from supersaturations lower than that needed for the mechanism described by the B+S model. The BCF model can be expressed mathematically as Eqs. (23)–(25):

$$G_{\mathit{BCF}}=\frac{{\tau }_{3}T}{{\tau }_4}{\left(ln\left(\sigma \right)\right)}^{2}tanh\left(\frac{{\tau }_4}{Tln\left(\sigma \right)}\right)$$E23

$${\tau }_3=\frac{{\mathrm{\Gamma }}^{\ast }{D}_{\mathit{\text{surf}}}{V}_m}{{x_s}^2}$$E24

$${\tau }_4=\frac{19V_m{\gamma }_{\mathit{sl}}}{2k{x}_s}$$E25

where ${\tau }_3$ and ${\tau }_4$ are lumped parameters and all the other parameters have the same meanings as before [6].

Now that we have presented the B+S and BCF models, we can appreciate the simplicity of the power-law model. In the introduction, we mentioned that Quilló and coworkers [6] used the power law, B+S, and BCF models, separately and in different combinations, to model the crystal growth kinetics of a proprietary API and thereby attempted to ascertain the underlying growth mechanism(s). During the combination tests, they also tested for dominance of one model over the other, sudden mechanism crossover, and simultaneous mechanisms. They found that the growth process could be best described if both the B+S and BCF models were simultaneously active, in an additive manner, with a smooth crossover from initial surface nucleation dominated growth to screw dislocation (spiral growth) dominated growth.

5. Crystallization from the amorphous state

One of the ways to improve the solubility of a drug is to prepare an amorphous solid from it. The molecules inside these amorphous solids lack the long-range ordered packing of their crystalline counterparts, giving the amorphous solids higher thermodynamic activity and generally enhanced solubility. Methods to enhance the solubility of drugs have become a hot topic over recent years, and will likely remain a topic of interest, especially considering that approximately 36% of the 698 drugs currently available as immediate-release oral preparations and 60% of the 28912 new APIs under development are classified as biopharmaceutics classification system (BCS) class II or IV drugs [30, 31]. The BCS attempts to classify a drug into one of four categories, depending on its aqueous solubility and membrane permeability [32]. BCS class II and IV drugs have poor aqueous solubility but good membrane permeability, and poor aqueous solubility and poor membrane permeability, respectively. In other words, the market share of poorly soluble drugs is likely to increase in the future.

Currently, there are a limited number of amorphous preparations on the market [33]. A specific cause for concern is that the high thermodynamic activity of an amorphous solid, responsible for its enhanced solubility, also makes it thermodynamically unstable and likely to convert back to a more stable, but less soluble, crystalline form. A mechanistic understanding of this conversion and the resulting glass-to-crystal (GC) growth is an important topic in pharmaceutical sciences.

GC growth rates cannot be readily explained by the thermodynamic driving force responsible for crystallization from solution, as described previously [34]. Initially, it was believed that crystallization from the amorphous state was dependent on the storage conditions, such as the storage temperature relative to the glass transition temperature $\left(T_g\right)$, and the bulk properties of the amorphous solid, like α-relaxation processes [34, 35, 36, 37]. However, it is now generally accepted that GC growth can take place at temperatures well below $T_g$ and follows two mechanisms, namely fast surface crystal growth and slower growth in the bulk (interior) [31].

Both the surface and bulk crystallization processes are believed to be controlled by the self-diffusion of the molecules in the amorphous solid. Using indomethacin as model API, both Wu and Yu [38] and Swallen and Ediger [39] found that the crystal growth rates in the interior of their amorphous solids were proportional to the self-diffusion coefficients for temperatures close to $T_g$. It is believed that surface crystallization is so much faster than bulk crystallization, because of the increased molecular mobility at the surface. This was corroborated by Zhu and coworkers [40] who measured the surface smoothing of indomethacin glasses and found the self-diffusion on the surface to be at least one million times faster than in the interior.

With the mechanism of surface crystallization understood, immobilization of the molecules at the surface makes for an appealing target to delay the crystallization process. Yu and coworkers [41] coated indomethacin glasses with gold (10 nm) and two polyelectrolytes. They found that even a single layer of polyelectrolyte was enough to inhibit the growth of exiting crystals and that the molecular mobility of molecules at the surface of an amorphous solid can be sufficiently suppressed by a coating only a few nanometers thick. Their results suggest that nanocoating is a promising technique to stabilize amorphous solids.

6. Conclusions

The crystallization of small-molecule APIs plays an important role in the pharmaceutical industry. Because the properties of these crystals, that is, crystal habit, size distribution, and polymorphic form, influence many downstream processes, ranging from quality control testing to formulation; it is desirable to be able to consistently produce crystals with specific properties. To achieve this kind of control over the crystallization process, a thorough understanding of the underlying mechanisms is required. In this chapter, we looked at the thermodynamic driving force behind crystallization, and how it can be rearranged into a supersaturation driving force. Methods to estimate the degree of supersaturation were discussed. The mechanisms behind crystal nucleation and growth, and the mathematical models describing these mechanisms, were discussed and examples of how these mechanisms can be used as intervention points to control the properties of the resulting crystals were given. We ended this chapter with a look at amorphous solids, which have a natural tendency to crystallize back to a more stable, but less soluble form, and saw that immobilization of the molecules at the surface of these solids with even a single layer of the polymer was enough to stabilize the amorphous solid. With an increasing movement in the pharmaceutical industry toward streamlining manufacturing processes through control, the techniques discussed in this chapter might see even more general use in the future.