Challenges in the use of atomistic simulations to predict solubilities of drug-like molecules

Abstract

Background: Solubility is a physical property of high importance to the pharmaceutical industry, the prediction of which for potential drugs has so far been a hard task. We attempted to predict the solubility of acetylsalicylic acid (ASA) by estimating the absolute chemical potentials of its most stable polymorph and of solutions with different concentrations of the drug molecule. Methods: Chemical potentials were estimated from all-atom molecular dynamics simulations. We used the Einstein molecule method (EMM) to predict the absolute chemical potential of the solid and solvation free energy calculations to predict the excess chemical potentials of the liquid-phase systems. Results: Reliable estimations of the chemical potentials for the solid and for a single ASA molecule using the EMM required an extremely large number of intermediate states for the free energy calculations, meaning that the calculations were extremely demanding computationally. Despite the computational cost, however, the computed value did not agree well with the experimental value, potentially due to limitations with the underlying energy model. Perhaps better values could be obtained with a better energy model; however, it seems likely computational cost may remain a limiting factor for use of this particular approach to solubility estimation. Conclusions: Solubility prediction of drug-like solids remains computationally challenging, and it appears that both the underlying energy model and the computational approach applied may need improvement before the approach is suitable for routine use.

Polarization Corrections and the Hydration Free Energy of Water

“…In particular, the solvation properties of water are of great interest for a vast range of applications. The key thermodynamic property to describe solvation processes in water is the Gibbs free energy of hydration ( Δ G Hyd ), and the calculation of Δ G Hyd of organic molecules is a huge area of research in computational chemistry. − In this paper, we set out to perform a systematic test of several classical nonpolarizable water models for their ability to predict the free energy of hydration of water, which is directly related to, but distinct from, its free energy of vaporization . In doing so, we uncovered what we believe to be a fundamental inconsistency in the way polarization effects are currently handled in the calculation of phase-change energies and free energies.…”

Prediction of the Relative Free Energies of Drug Polymorphs above Zero Kelvin

“…The free energy difference between a fully interacting crystal and the EC state is derived by molecular simulation using a series of overlapping intermediate states along a thermodynamic pathway (Figure b). In previous ECMs, the simulations involved a two-step process including adding harmonic restraints to all atoms in a single step followed by removing all interactions. ,, In crystals of small rigid molecules, this process results in free energies that converge relatively quickly. However, adding all-atom restraints to the entire crystal in a single step can be problematic for more complex systems relevant to the pharmaceutical industry, where disordered atoms may undergo large fluctuations away from their preferred lattice positions, and on slow transition time scales .…”

“…For arbitrary temperatures, one could run separate ECM calculations for each temperature of interest. However, an alternative approach is to run a series of simulations for each crystal over a range of temperatures in the NPT ensemble and compute their reduced free energy using standard techniques like Thermodynamic Integration (TI) or the Multistate Bennett Acceptance Ratio (MBAR). , The dimensionalized relative Gibbs free energy difference is then recovered from the reduced free energies according to eq . where Δ f ij ( T ) is the reduced free energy difference between crystal i and j at temperature T and Δ G ij ( T ref ) is the Gibbs free energy difference at a reference temperature, T ref , computed from a single ECM calculation. All free energy calculations in this work, including the reduced free energies along temperatures and alchemical free energy changes along the thermodynamic path, are computed using MBAR.…”

“…In previous thermodynamic stability studies of polymorphic systems, it was observed that discrepancies between experiment and theory with point-charge force fields derive primarily from enthalpic rather than entropic differences. Furthermore, enthalpy differences between polymorphs tend to be nearly constant across temperature ranges . Therefore, after computing the free energy difference from eq , we shift the free energy profiles using a constant correction term (per polymorph) based on DFT lattice energies in conjunction with a harmonically derived zero point vibrational energy (ZPVE) as given by eq . where Δ E ij DFT ( T ) and Δ E ij MM are the DFT and MM minimized energy difference between the lattice energies of forms i and j , Δ E ij ZPVE is the harmonically derived ZPVE correction at zero Kelvin, and ΔΔ G ij MM ( T ) is the relative free energy change between forms i and j .…”

“…Then by evaluating the free energy of the nonphysical Einstein Crystal (EC) state, the free energy difference of the two physical polymorphs can be computed. Several variants of this method have been proposed in recent years, with subtle differences in the definition of the analytical state and the precise thermodynamic path traversed. − Here, we adopt a methodology similar to the pseudo supercritical path (PSCP) method, which was originally designed to compute the melting temperature of organic crystals and was later extended to rank the finite-temperature stability of polymorph pairs . The advantage of this method is that it is possible to separate the enthalpic and entropic contributions to the free energy, thereby making it possible to combine the accuracy of modern electronic structure methods such as DFT with the ability to effectively sample an ensemble of configurations using classical molecular dynamics.…”