simulation of dissolution

Noyes–Whitney Dissolution Engine

Biopharmaceutics & Drug Release Kinetics

Governing Transport ODE Sink Conditions Active (< 15% Cs)

$$\frac{dC}{dt} = \frac{D \cdot A(t)}{h \cdot V} \cdot (C_s - C(t))$$

Formulation Presets:

D: Diffusion Coeff

Solute mobility across boundary layer

A(t): Surface Area

Solute-solvent interfacial contact

h: Stagnant Layer

Diffusion thickness (modulated by stirring)

Cs: Saturation Sol.

Thermodynamic equilibrium solubility

C(t): Bulk Conc.

Instantaneous dissolved concentration

Cs - C: Driving Force

Chemical potential gradient driving dissolution

Physical Parameters

SI & Bio units

Saturation Solubility ($C_s$) 20.0 µg/mL

Poorly Soluble (2 µg/mL) Moderate (100 µg/mL)

Initial Particle Radius ($r_0$) 5.0 µm

Sub-micron (0.5 µm) Coarse Granule (80 µm)

Boundary Layer Thickness ($h$) 25.0 µm

Turbulent / High RPM (5 µm) Gentle / Static (60 µm)

Diffusion Coefficient ($D \times 10^{-6}$) 5.5 cm²/s

Typical small molecule aqueous diffusivity is 4–8 × 10⁻⁶ cm²/s

Drug Dose ($M_0$)

mg

Fluid Volume ($V$)

mL

Simulation Horizon 60 min

Overlay Constant Surface Area profile Show Saturation Solubility ($C_s$) line

Particle & Stagnant Layer ($h$)

r(t): 5.0 µm

Inner: Solid API • Outer ring: Diffusion Layer $h$

Drug concentration declines smoothly from $C_s$ at the particle interface across stagnant boundary thickness $h$ into bulk solution $C(t)$.

Dissolution at Horizon 0.0% 0.00 mg dissolved

Initial Rate ($dC/dt_0$) 0.00 µg / (mL · min)

Time to 80% ($T_{80}$) -- Pharmacopeial threshold

Specific Surface Area 0.00 m² / g API

Drug Concentration Profile $C(t)$

Shrinking Particle (Spherical ODE) vs. Idealized Constant Surface Area

Shrinking Particle (Realistic) Constant Area

Particle Radius $r(t)$ & Area Decay

Micrometers (µm)

Instantaneous Rate $dC/dt(t)$

µg / (mL · min)

Theoretical Foundations & Biopharmaceutics Formulation Levers

How physicochemical parameters dictate in vitro dissolution testing and oral bioavailability (BCS Class II drugs).

1. Surface Area ($A$) & Micronization

The dissolution rate is directly proportional to total interfacial contact area $A$. For spherical particles of radius $r_0$ and density $\rho$, specific surface area scales as: SSA = 3 / (ρ · r₀) Reducing particle radius by a factor of 10 increases the available surface area tenfold for the same mass dose, dramatically accelerating the initial dissolution velocity. This is the biopharmaceutical rationale for micronization, wet bead milling, and nanosuspensions.

2. Hydrodynamics & Boundary Layer ($h$)

Surrounding every dissolving particle is a microscopic stagnant fluid film across which solute molecules must diffuse purely by Brownian motion. According to Levich hydrodynamic theory, agitation speed ($\omega$) compresses this layer: h ∝ ω^(-1/2) Higher gastrointestinal motility or dissolution paddle speeds (e.g. 50 vs. 100 RPM in USP Apparatus II) thin this barrier, elevating $dC/dt$.

3. Driving Force $(C_s - C)$ & Sink Conditions

Dissolution rate is driven by chemical chemical disequilibrium. As solute accumulates in the bulk volume, $C(t) \to C_s$, reducing driving force $(C_s - C)$ towards zero. In biopharmaceutics, sink conditions are defined when $C(t) < 0.1 \sim 0.2 \, C_s$, preserving maximal driving force. Strategies like surfactant micellization, salt formation, or amorphous solid dispersions (ASDs) generate temporary "spring and parachute" supersaturation by elevating $C_s$.

Coupled ODE Formulation: Shrinking Sphere Kinetics

A common error in naive dissolution models is treating surface area $A$ as a static constant. As mass dissolves, the particle shrinks. Assuming $N$ identical spherical particles of density $\rho$, the single particle mass is $m_p = \frac{4}{3}\pi r^3 \rho$ and surface area is $a_p = 4\pi r^2$. Taking the time derivative:

Particle Radius Shrinkage ODE: $$\frac{dr}{dt} = - \frac{D}{\rho \cdot h} \cdot \left(C_s - C(t)\right) \quad (r > 0)$$

Bulk Concentration ODE: $$\frac{dC}{dt} = \frac{N \cdot 4\pi r(t)^2 \cdot D}{h \cdot V} \cdot \left(C_s - C(t)\right)$$

When $C \ll C_s$ (steady sink condition), integrating the radius equation directly yields the classic Hixson–Crowell Cube-Root Law: $M_0^{1/3} - M(t)^{1/3} = \kappa \cdot t$.

Standalone Python Implementation (SciPy & Matplotlib)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# =====================================================================
# NOYES-WHITNEY DRUG DISSOLUTION SIMULATION (COUPLED ODE)
# =====================================================================

# 1. Parameter Specifications for a Poorly-Soluble Drug (e.g. Griseofulvin)
Dose_mg = 10.0            # Total drug dose (mg)
V_mL = 500.0              # Dissolution medium volume (mL)
Cs_ug_mL = 20.0           # Saturation solubility (µg/mL)
r0_um = 5.0               # Initial particle radius (µm)
h_um = 25.0               # Diffusion boundary layer thickness (µm)
D_cm2_s = 5.5e-6          # Diffusion coefficient (cm²/s)
rho_g_cm3 = 1.25          # Solid drug density (g/cm³)
t_span_min = 60.0         # Simulation duration (minutes)

# Unit Conversions into SI / consistent CGS units
Cs = Cs_ug_mL * 1e-3      # mg/mL
r0 = r0_um * 1e-4         # cm
h = h_um * 1e-4           # cm
D = D_cm2_s * 60.0        # cm²/min (converted to minutes)
rho = rho_g_cm3 * 1e3     # mg/cm³

# Calculate number of particles (N) based on initial dose and spherical geometry
single_particle_volume_0 = (4.0 / 3.0) * np.pi * (r0 ** 3)
single_particle_mass_0 = single_particle_volume_0 * rho  # mg
N_particles = Dose_mg / single_particle_mass_0
initial_total_area = N_particles * (4.0 * np.pi * (r0 ** 2)) # cm²

# ---------------------------------------------------------------------
# Scenario 1: Realistic Coupled ODE (Shrinking Spherical Particles)
# State vector y = [r, C] where r is particle radius (cm), C is conc (mg/mL)
# ---------------------------------------------------------------------
def shrinking_particle_odes(t, y):
    r, C = y
    
    # Boundary condition: if particle is completely dissolved or solution saturated
    if r <= 0.0 or C >= Cs:
        return [0.0, 0.0]
    
    driving_force = max(0.0, Cs - C)
    
    # Radius shrinkage rate: dr/dt = - (D / (rho * h)) * (Cs - C)
    dr_dt = - (D / (rho * h)) * driving_force
    
    # Mass balance for bulk concentration:
    # dC/dt = (N * 4 * pi * r^2 * D / (h * V)) * (Cs - C)
    current_area = N_particles * (4.0 * np.pi * (r ** 2))
    dC_dt = (D * current_area / (h * V_mL)) * driving_force
    
    return [dr_dt, dC_dt]

# ---------------------------------------------------------------------
# Scenario 2: Constant Surface Area Approximation (Static A0)
# ---------------------------------------------------------------------
def constant_area_ode(t, y):
    C = y[0]
    if C >= Cs:
        return [0.0]
    driving_force = max(0.0, Cs - C)
    dC_dt = (D * initial_total_area / (h * V_mL)) * driving_force
    return [dC_dt]

# Time evaluation points
t_eval = np.linspace(0, t_span_min, 500)

# Solve coupled system
sol_shrinking = solve_ivp(shrinking_particle_odes, [0, t_span_min], [r0, 0.0], t_eval=t_eval)
sol_constant = solve_ivp(constant_area_ode, [0, t_span_min], [0.0], t_eval=t_eval)

# Extract solutions (convert concentration back to µg/mL)
C_shrinking_ug_mL = sol_shrinking.y[1] * 1e3
r_shrinking_um = sol_shrinking.y[0] * 1e4
C_constant_ug_mL = sol_constant.y[0] * 1e3

# ---------------------------------------------------------------------
# Plotting Results
# ---------------------------------------------------------------------
plt.figure(figsize=(10, 5), dpi=120)
plt.plot(t_eval, C_shrinking_ug_mL, 'b-', lw=2.5, label="Shrinking Particle Model (Realistic)")
plt.plot(t_eval, C_constant_ug_mL, 'r--', lw=1.8, label="Constant Surface Area ($A_0$)")
plt.axhline(Cs_ug_mL, color='gray', linestyle=':', label=f"Saturation Solubility $C_s$ ({Cs_ug_mL} µg/mL)")
plt.title("Noyes–Whitney Dissolution Kinetics: Shrinking vs Constant Area", fontsize=12)
plt.xlabel("Time (minutes)", fontsize=10)
plt.ylabel("Bulk Drug Concentration (µg/mL)", fontsize=10)
plt.grid(True, linestyle='--', alpha=0.6)
plt.legend(loc='lower right')
plt.tight_layout()
plt.show()