Coupling Nucleation and Growth in Precipitation Models

Why Nucleation and Growth Must Be Coupled

Supersaturation can produce solid drug by two linked processes: new crystals appear, and crystals already present grow. A model that includes growth but omits new nuclei may underpredict the number of small particles formed. A model that includes nucleation but treats the crystals as having no size-dependent mass cannot calculate how much dissolved drug their growth consumes. Coupling the processes makes the changing particle population and its effect on dissolved concentration explicit.

The model below describes precipitation in a well-mixed liquid of fixed volume. It tracks spherical crystals by radius and uses a precipitation-only scope: crystals do not dissolve in this model. That assumption is most suitable while the solution remains supersaturated. The model does not specify how a changing supersaturation might later cause redissolution.

Definition: Let $C(t)$ be dissolved-drug concentration and $C_s$ its equilibrium solubility in the modeled liquid phase. The supersaturation ratio is $S(t)=C(t)/C_s$. Let $N(r,t)$ be the number of crystals per unit liquid volume per unit radius at radius $r$ and time $t$; $J(S)$ is the rate of formation of nuclei per unit liquid volume; and $G(r,S)=dr/dt$ is the radial crystal-growth rate. Thus, $N(r,t)\,dr$ is the number of crystals per unit liquid volume with radii between $r$ and $r+dr$.

A useful illustrative pair of rate laws is a power law in the supersaturation excess:

J(S)=k_n [S-1]_+^b, \qquad G(r,S)=k_g [S-1]_+^g, \qquad [x]_+=\max(x,0).

Here $k_n$ has units of nuclei per unit liquid volume per unit time, $k_g$ has units of radius per unit time, and $b$ and $g$ are dimensionless exponents. The displayed growth law has no explicit size dependence, but a more detailed model can use $G(r,S)$. The power laws are examples of kinetic laws, not universal mechanisms. In either case, this precipitation-only formulation sets nucleation and growth to zero for $S\leq 1$; it does not represent dissolution there.

Representing the Evolving Crystal Population

Suppose nuclei enter the growing population at a specified radius $r_0>0$. For radii above $r_0$, crystals move through size space at velocity $G$. If crystals are neither lost nor broken, their population satisfies the size-structured balance

\frac{\partial N(r,t)}{\partial t} + \frac{\partial}{\partial r}\left(G(r,S(t))N(r,t)\right)=0, \qquad r>r_0.

New nuclei enter through the lower-radius boundary. When growth is positive, the influx condition is $G(r_0,S)N(r_0,t)=J(S)$. The two sides have units of number per unit liquid volume per unit time. This condition specifies a particle flux into the size domain, not a concentration of dissolved drug.

The size distribution is needed because the mass of a spherical crystal depends on its radius. If $\rho_s$ is the solid density, the mass of one crystal of radius $r$ is $m(r)=4\pi\rho_s r^3/3$. The total suspended solid mass per unit liquid volume is therefore

M_s(t)=\int_{r_0}^{\infty}m(r)N(r,t)\,dr.

These equations assume a well-mixed liquid, a common solubility, spherical crystals of fixed density, and precipitation as the only process changing dissolved drug. They also assume the population has no aggregation, fragmentation, settling, or outflow. Such processes require additional terms; they cannot be represented by changing $J$ or $G$ alone.

How the Population Depletes Dissolved Drug

A nucleus added at the boundary already has mass $m(r_0)$. Meanwhile, every crystal in the population gains mass as its radius increases. Consequently, the precipitation rate contains both a birth contribution and a growth contribution. For fixed liquid volume and no other source or sink of drug, the dissolved concentration obeys

\frac{dC}{dt} = -\left[ m(r_0)J(S) + \int_{r_0}^{\infty} \frac{dm(r)}{dr}G(r,S)N(r,t)\,dr \right].

The first term is the mass carried by newly formed nuclei per unit liquid volume per unit time. The integral is the mass gained per unit time by the existing population. Since $dm/dr=4\pi\rho_s r^2$, the growth contribution weights each crystal by its surface area, its radial growth rate, and its abundance at that size. The minus sign indicates depletion of dissolved drug.

Derivation: Multiply the population balance by $m(r)$ and integrate over $r\geq r_0$. Integration by parts, with no particle flux at infinite radius, gives $$ \frac{dM_s}{dt} = m(r_0)G(r_0,S)N(r_0,t) + \int_{r_0}^{\infty}\frac{dm}{dr}G(r,S)N(r,t)\,dr. $$ Using the boundary influx $G(r_0,S)N(r_0,t)=J(S)$ yields the two terms in the solid-mass rate. If precipitation transfers drug only from dissolved to solid form, then $dC/dt=-dM_s/dt$, which gives the stated concentration equation.

This derivation shows why multiplying $J$ by a growth rate, or using only a bulk growth term, does not generally give the correct depletion rate. Nucleation supplies particles at a defined initial mass; subsequent growth acts on all crystals, and its contribution depends on the current size distribution.

Worked Example: Evaluating Both Kinetic Rates

Consider hypothetical kinetic parameters $k_n=2.00\times10^8$ nuclei mL$^{-1}$ min$^{-1}$, $b=2$, $k_g=0.400$ μm min$^{-1}$, and $g=1.5$. At $S=1.50$, calculate the nucleation and growth rates.

The supersaturation excess is $S-1=0.50$. Substituting into the nucleation law gives

J=(2.00\times10^8\ \text{nuclei mL}^{-1}\text{ min}^{-1})(0.50)^2 =5.00\times10^7\ \text{nuclei mL}^{-1}\text{ min}^{-1}.

The growth rate is

G=(0.400\ \text{μm min}^{-1})(0.50)^{1.5} =(0.400)(0.3536) =0.1414\ \text{μm min}^{-1}.

Thus, under these hypothetical conditions, nuclei are being added while crystals increase in radius by about $0.1414$ μm per minute. Both rates are positive because $S>1$. These rate calculations alone do not determine the concentration decline: the current crystal-size distribution and the mass of a new nucleus are also needed.

Worked Example: Calculating Joint Dissolved-Drug Depletion

Use the rates from the preceding example. Suppose, hypothetically, the solid density is $1.20\times10^{-9}$ mg μm$^{-3}$, the nucleus radius is $r_0=0.0500$ μm, and the existing population is monodisperse: $1.00\times10^8$ crystals mL$^{-1}$ at radius $1.00$ μm. Assume the same $S=1.50$ applies to all crystals. Calculate the instantaneous depletion rate.

First, the mass of a nucleus at the boundary is

m(r_0)=\frac{4\pi}{3}(1.20\times10^{-9})(0.0500)^3 =6.283\times10^{-13}\ \text{mg per crystal}.

The nucleation contribution to precipitation is therefore

m(r_0)J =(6.283\times10^{-13})(5.00\times10^7) =3.142\times10^{-5}\ \text{mg mL}^{-1}\text{ min}^{-1}.

For one existing crystal at $r=1.00$ μm, the growth-related mass gain is

\frac{dm}{dr}G =4\pi(1.20\times10^{-9})(1.00)^2(0.1414) =2.132\times10^{-9}\ \text{mg min}^{-1}.

Multiplying by the existing number concentration gives a growth contribution of $0.2132$ mg mL$^{-1}$ min$^{-1}$. Hence

\frac{dC}{dt} =-(0.2132+0.00003142) \approx -0.2133\ \text{mg mL}^{-1}\text{ min}^{-1}.

Growth of the existing crystals dominates this particular rate calculation; nucleation still contributes, but the mass of each new nucleus is small. This is an instantaneous rate, not a prediction that concentration will continue falling linearly at this rate. As $C$, $S$, and the population change, the rates generally change as well.

New Nuclei Become a Distribution of Sizes

A continuous nucleation rate creates crystals of different ages at any observation time. Even if all nuclei enter at the same radius, older nuclei have had more time to grow. When growth depends on time through supersaturation, a particle born at time $\tau$ follows the characteristic equation $dr/dt=G(r,S(t))$, starting at $r(\tau)=r_0$. Its current radius is therefore determined by the history of $S$ since birth, not just the current supersaturation.

Worked Example: Following a Newly Formed Cohort

For a hypothetical illustration of population evolution, suppose $S$ is held at 1.50 for two minutes, so that $J=5.00\times10^7$ nuclei mL$^{-1}$ min$^{-1}$ and $G=0.1414$ μm min$^{-1}$ remain constant. Nuclei enter at $r_0=0.0500$ μm. Find the number and range of radii of nuclei formed during this interval, and their mean radius at the end.

The number added per millilitre over two minutes is

(5.00\times10^7\ \text{nuclei mL}^{-1}\text{ min}^{-1})(2.00\ \text{min}) =1.00\times10^8\ \text{nuclei mL}^{-1}.

A nucleus formed at time $\tau$ has age $2-\tau$ minutes at the end of the interval. With constant growth, its radius is $r_0+G(2-\tau)$. Nuclei formed at the end remain at $0.0500$ μm; those formed at the start reach

0.0500+(0.1414)(2.00)=0.3328\ \text{μm}.

Births occur uniformly over the interval because $J$ is constant. Their ages at the endpoint range uniformly from 0 to 2 minutes, so the mean age is $1.00$ minute. Their mean radius is therefore $0.0500+(0.1414)(1.00)=0.1914$ μm. The newly formed population spans radii from $0.0500$ to $0.3328$ μm; representing all of it as particles at the mean radius would discard that spread. This example holds $S$ fixed to make the population calculation transparent; it is not a closed-liquid concentration prediction.

Model Choices, Diagnostics, and Limitations

A size-resolved population balance is one option for coupling the processes. A simpler alternative is a lumped precipitation law that directly specifies the rate at which dissolved drug becomes solid. Such a law can be useful when only dissolved concentration is of interest and particle size has little influence on the application. However, it does not separately predict how many crystals form, how they grow, or how their sizes affect later growth-related depletion. A representative-radius model retains some size information with fewer state variables, but can misrepresent a broad or multimodal population. The population balance is more informative when size distribution matters, at the cost of additional numerical and parameter requirements.

Several checks help identify a model or implementation that is going wrong:

Check the supersaturation regime. Under the stated laws, $S\leq1$ gives $J=G=0$. Positive nucleation or growth in that regime indicates a coding or model-definition mismatch.
Check units. $m(r_0)J$ and the growth integral must both have units of mass per liquid volume per time. The boundary flux $GN$ must have the same units as $J$.
Check population movement. If $G>0$, crystals should move toward larger radii. A negative population density or unexplained loss of particles suggests a discretization or boundary-flux problem.
Check concentration response. With positive nucleation and growth and no other processes, $dC/dt\leq0$. A computed increase in dissolved concentration is inconsistent with this precipitation-only model.
Check numerical resolution. If newly formed nuclei cross several size bins in one time step, or supersaturation changes sharply within a step, a coarse grid or time step may distort the boundary influx and depletion rate. Refine the grid or step and check whether predictions stabilize.

The model also has substantive limitations. Spherical shape, constant density, and a single well-mixed dissolved concentration may be poor approximations for some formulations. The kinetic laws may not transfer across media or conditions without evidence. A population balance does not by itself establish that the kinetic parameters can be uniquely inferred from concentration data. Most importantly, the precipitation-only scope is inadequate once redissolution or other dissolution processes become relevant; those processes need to be represented explicitly rather than interpreted as negative nucleation.

Key takeaway: Nucleation determines the influx of crystals at the smallest modeled radius; growth moves the entire population through size space. Dissolved drug is depleted both to supply the mass of new nuclei and to supply the mass gained as existing crystals grow.

Check Your Understanding

Answer each question, then compare with the solutions below.

If $S=1$ under the power-law model, what are $J$ and $G$, and what does that imply for precipitation in this formulation?
Why does the dissolved-drug depletion equation include a term $m(r_0)J$ as well as a growth integral?
Suppose a crystal population has $N=2.00\times10^7$ crystals mL$^{-1}$ at $r=2.00$ μm, with $\rho_s=1.00\times10^{-9}$ mg μm$^{-3}$ and $G=0.100$ μm min$^{-1}$. What is its growth-related depletion rate?
If a constant nucleation rate operates for three minutes, why are the new crystals generally not all the same radius at the end?

Solutions. (1) Since $[S-1]_+=0$, both rates are zero. This model predicts no new nuclei or crystal growth at $S=1$; it does not describe dissolution. (2) The boundary influx introduces nuclei with finite radius and hence nonzero initial mass, giving $m(r_0)J$. The integral accounts for mass subsequently gained by all crystals as their radii increase. (3) Per crystal, $4\pi\rho_s r^2G=4\pi(1.00\times10^{-9})(2.00)^2(0.100)=5.027\times10^{-9}$ mg min$^{-1}$. Multiplying by $2.00\times10^7$ crystals mL$^{-1}$ gives $0.1005$ mg mL$^{-1}$ min$^{-1}$. (4) They have different ages: crystals formed earlier have had more time to grow. Even with constant $G$, the earliest and latest births differ in growth time by three minutes, so their radii differ by $3G$.