Predicting Aqueous Solubility of Pharmaceutical Cocrystals Using Constant Chemical Potential Enhanced Sampling Simulations

Abstract We present an integrated molecular simulation framework that combines molecular dynamics, enhanced sampling, and constant chemical potential simulations to uncover growth and dissolution mechanisms of molecular cocrystals and predict equilibrium solubility, using the theophylline–aspirin (THE–ASP) cocrystal as a prototypical system. Solubility is determined by sampling the growth and dissolution of ASP and THE molecules at a dimeric kink site on the crystal surface and computing the free-energy differences between their dissolved and crystalline states as a function of solution composition. These free-energy differences are interpreted as effective chemical potentials of the respective molecular species, and the solubility is identified as the solution composition at which the sum of the effective chemical potentials of ASP and THE equals zero. To address the large conformational and orientational freedom inherent to these molecules, we introduce collective variables that incorporate atomic density, molecular orientation, and solvent coordination, which, when used in Well-tempered Metadynamics simulations, enable efficient sampling of growth and dissolution processes and the calculation of well-converged free-energy surfaces. The calculated free energy profiles reveal asymmetric thermodynamic driving forces, with aspirin favoring incorporation into the crystal and theophylline preferentially stabilizing in solution, providing molecular-level insight into cocrystal growth dynamics. The predicted equilibrium solubility of the THE–ASP cocrystal, χ* = 0.00063 ± 0.00014, lies between the reported aqueous solubilities of the individual components (ASP ≈ (3–5) × 10–4; THE ≈ (6–8) × 10–4 at 25 °C).