(M1230-06-39) A Breakthrough Predictive Model for Rapid and Accurate Solubility Parameter Prediction of Lipids for Lipid-Based Drug Delivery Systems

Purpose: Preformulation and formulation epitomize the phases involved in pharmaceutical product development, and determining the solubility is the cornerstone of this process. The Hildebrand solubility Parameter (δ)(1)(Figure1:Eq.1), a quantitative metric of cohesive energy density (MPa^½), delineates a substance’s dissolution tendency via the dictum “like dissolves like”. This parameter integrates electrostatic interactions (δElec), and van der Waals forces (δVdW), offering profound insights into drug-excipient miscibility and solubility. Despite its utility, traditional experimental methods are laborious, prompting the adoption of molecular dynamics (MD) simulations for computational estimation. However, MD simulations grapple with numerous challenges: high computational demands (~25 days on a CPU or ~2 days on a GPU, per simulation), computational architecture, memory constraints, and intricate force field requirements, rendering them less pragmatic for large-scale applications. To circumvent these issues, a mathematical predictive model is proposed for lipid excipients, leveraging physicochemical descriptors such as δElec, δVdW, and simulated specific heat capacity. This model would eliminate the need for computationally onerous MD simulations, enabling rapid but accurate predictions of solubility parameters for lipids and other pharmaceutical excipients. This approach streamlines the rational design and optimization of lipid-based delivery systems—such as lipid nanoparticles and self-emulsifying drug delivery systems (SEDDS), which are crucial for improving the absorption of poorly soluble drugs.

Methods: MD simulations were used to predict key physicochemical properties of the lipids (Fig.2) using Schrödinger Materials Maestro(2) (Schrödinger Inc., NY, USA) Desmond(3) engine, enabling prediction of bulk properties of the lipids. Structured data files (SDF) for each lipid were subjected to Low-mode sampling (LMOD) was employed for conformational exploration, leveraging eigenvectors to identify flexible degrees of freedom. Parameters included 1000 steps per molecule, 100 steps per rotatable bond, and an energy window of 50 kJ/mol. Structures differing by >0.5 Å in atomic positions were retained for obtaining the global minimum energy conformation, which was selected for further simulations. A single-component amorphous disordered system of 400 lipid molecules was constructed in vacuum using a constant dielectric model (ε=1.0), under periodic boundary conditions (PBC), to isolate lipid interactions. Nonbonded interactions were modelled using cutoff distances: 7.0Å for van der Waals, 12.0Å for electrostatic, and 4.0Å for hydrogen bonding. The system was equilibrated using a six-step relaxation protocol, followed by a 100 ns NPT-ensemble MD simulation at 298K and 1.01325bar, with a 2fs timestep. A Nosé-Hoover thermostat and Martyna-Tobias-Klein barostat ensured temperature and pressure control. Post-equilibration, bulk properties such as density, cohesive energy, heat of vaporization, and solubility parameters were extracted from the final 20% of the trajectory. This protocol was applied to 76 lipids, divided into training and test sets using RandomLists. Principal Component Analysis and stepwise regression identified significant predictors, which were analyzed statistically.

Results: Principal component analysis (PCA) distilled the critical physicochemical descriptors to Electrostatic forces (δElec), Van der Waals forces (δVdW), and Simulated specific heat capacity, which were subsequently harnessed for computing solubility parameters (δ). Stepwise multiple regression analysis yielded several models, among which the optimal model exhibiting the highest adjusted R²=98.21% (adj. R²=0.9821), was selected for further investigation. The derived predictive equation demonstrated robustness in estimating solubility parameters (Fig1:Eq.2). The disparity between experimental and predicted datasets exhibited a mean error of 4.3895%, indicative of satisfactory congruence. The Shapiro-Wilk test for normality revealed significant normality in both training (W=0.93105; p-value=0.05855) and test datasets (W=0.96592; p-value=0.455), facilitating parametric assumptions. Subsequent application of the parametric two-sample t-test (α=0.05) yielded p-value=0.1365, suggesting no statistically significant difference between the two groups. Furthermore, Pearson’s product-moment correlation coefficient (r=0.879; 95%CI=0.757-0.942)(Fig.3a) underscored a strong monotonic relationship between the datasets, reinforcing their consistency. The R²=77.3% (adj. R²=76.5%)(Fig.3b) and Root Mean Square Deviation (RMSD) value computed to quantify similarities was 1.271, corroborating the model’s strong predictive accuracy. Overall, these values represent exceptional performance for a predictive model, significantly accelerating preformulation studies and improving the accuracy of solubility predictions, particularly for poorly water-soluble drugs.

Conclusion: The development and implementation of this mathematical predictive model represent a significant frog-leap in pharmaceutical research, poised to redefine the landscape of preformulation and formulation studies. By seamlessly integrating physicochemical descriptors such as δElec, δVdW, and simulated specific heat capacity, this model offers unparalleled precision in predicting solubility parameters for lipid excipients while entirely circumventing the computational bottlenecks associated with traditional MD simulations. This transformative approach not only expedites the rational design and optimization of lipid-based drug delivery systems—such as lipid nanoparticles and SEDDS—but also enhances their efficacy in addressing the challenges posed by poorly water-soluble drugs. This research highlights a transformative predictive model that delivers rapid, accurate, and cost-effective solubility predictions, thus revolutionizing pharmaceutical development by enhancing efficiency, scalability, and accessibility while maintaining scientific rigor and economic feasibility.

References: (1) J. H. Hildebrand, Solubility of Non-Electrolytes, 2nd ed. (Reinhold, New York, 1936; Russian translation: GONTI, Moscow, 1938).
(2) Schrödinger Materials Maestro-Materials Science; New York, NY, 2019.
(3) Schrödinger Maestro-Desmond Interoperability Tools; New York, NY, 2019.

Acknowledgements: We sincerely thank Dr. Sriramakamal Jonnalagadda, Dr. Jasmin Monpara, and Dr. Rijo John for their invaluable guidance and support. We also extend our gratitude to all employees of the Industrial Pharmacy Laboratory (IPhL) for their dedication and contributions to this research.


Fig 1: (Eq.1) Hildebrand Solubility Parameter, where, Hvap: Heat of Vaporization, R: Universal Gas Constant, T: Temperature (in Kelvin), and Vm: molar volume of components.
(Eq.2) Equation derived from the Predictive Model, where, :Hildbrand Solubility Parameter, VdW:Van der Waals forces, Elec:Electrostatic forces, Hspec:Simulated specific Heat Capacity.


Fig 2: Trajectory snapshots of the amorphous disordered system after the six-step relaxation and equilibration for Cholesterol (a part of the training set) at (a) 0 nanoseconds; (b) 50 nanoseconds; and (c) 100 nanoseconds of molecular dynamics. Similarly, trajectory snapshots of the amorphous disordered system after the six-step relaxation and equilibration for Dipalmitoyl phosphatidylserine (DPPS) (a part of the test set) at (d) 0 nanoseconds; (e) 50 nanoseconds; and (f) 100 nanoseconds of molecular dynamics.


Figure 3: (A) Matrix plot of the Pearson Product-Moment Correlation, yielding a linear positive correlation (r) of 0.879 (95% Confidence Interval = 0.757 - 0.942), demonstrating a substantial correlation between experimental and predictive model results. (B) Scatterplot of the relation between the simulated and predicted values of solubility parameters, yielding an adj. R2 of 76.5% (0.765), substantial for predictive models.