Computation of drug solvation free energy in supercritical CO${}_2$: Alternatives to all-atom computer simulations

Abstract

Despite the modern level of development of computational chemistry methods and technological progress, fast and accurate determination of solvation free energy remains a huge problem for physical chemists. In this paper, we describe two computational schemes that can potentially solve this problem. We consider systems of poorly soluble drug compounds in supercritical carbon dioxide. Considering that the biggest contribution among all intermolecular interactions is made by van der Waals interactions, we model solute and solvent particles as coarse-grained ones interacting via the effective Lennard-Jones potential. The first proposed approach is based on the classical density functional theory and the second one relies on molecular dynamics simulation of the Lennard-Jones fluid. Sacrificing the precision of the molecular structure description while capturing the phase behavior of the fluid with sufficient accuracy, we propose computationally advantageous paths to obtaining the solvation free energy values with the accuracy satisfactory for engineering applications. The agreement reached between the results of such coarse-graining models and the experimental data indicates that the use of the all-atom molecular dynamic simulations for the studied systems seems to be excessive.

Introduction

Solvation free energy is an extremely important thermodynamic quantity for a wide range of physico-chemical applications, including calculation of solubility values, estimation of partition coefficients, infinite dilution activity coefficients and other parameters. Nowadays there are a number of computational methods developed to determine solvation free energy values. Regarding the systems that include poorly soluble drugs and drug-like compounds it is rather popular among scientists to rely on the results of the all-atom molecular dynamics (AAMD) simulations based on alchemical free energy methods [1], [2], [3], [4], [5], [6]. Such approaches model a set of nonphysical intermediate states to compute the values of free energy of solute molecule transfer from the gas phase to the solution. More details regarding different alchemical methods can be found elsewhere [7]. Despite their evident efficiency and similarly clear time expenditures and computational cost, these methods can have other drawbacks. They include high sensitivity to the force field choice and uncertainty arising from the choice of the partial atomic charges computation method. Both of these aspects of the AAMD calculation can lead to the inconsistency of the final solvation free energy values up to several kcal/mol [1], [8], [9]. At the same time, the procedure of the force-field parametrization is itself a nontrivial task and it is not always easy to justify some of the approaches physically. Thus, if the aim is to determine the free energy quickly, the utilization of the AAMD techniques is rather questionable. On the other hand, as it has been recently shown [10], [11], [12], [13], the results of the coarse-graining models using united-atom force fields are comparable in accuracy with those of the all-atom ones and in some cases even outperform them. Another approach to the estimation of drug solvation free energy in supercritical CO${}_2$ (scCO${}_2$), popular among chemical engineers, is based on applying cubic equations of state, for instance, the Peng-Robinson equation of state (PR EOS) for binary mixtures, which is a generalization of the classical van der Waals equation of state [14], [15]. However, the PR EOS contains several adjustable parameters, the determination of which requires an experimental data set (solubility isotherms, for instance). That is why the PR EOS is usually used to approximate the existing experimental data. Thus, in order to apply the PR EOS for solvation free energy estimation, it is necessary to perform preliminary experimental measurements that make this approach much more complicated. The next step toward a more robust description of the phase equilibrium is a family of SAFT (Statistical Associating Fluid Theory) type models [16]. As they include the contributions of different interaction terms and the parameters of these models hold a strict physical meaning, they have been successfully utilized to correlate and predict the solvation properties of a variety of solutes [17], [18], [19], [20], [21], [22]. Although there is a rather substantial decrease in the number of binary interaction parameters, when multicomponent systems are considered, as compared to the cubic EOS with advanced mixing rules, the whole parametrization procedure is still cumbersome and the choice of the best one is itself a nontrivial task [23]. High accuracy of the calculated solvation free energy values with an error less than 1 kcal/mol for sufficiently large data sets can be also achieved by using quantum mechanical calculations on continuum solvents (COSMO-based models, SMD) [24], [25], [26], [27], [28]. The cost of such accuracy, however, is the complex parametrization routine and the necessity of its replication for new systems, if they are dramatically different from the original learning set [29].

Quite promising is the method of solvation free energy estimation based on the classical density functional theory (cDFT). In papers [30], [31], the authors proposed a method for calculating hydration free energy of a set of hydrophobic solutes. This method can be considered as a coarse-grained one, since within its framework real molecules of both the solute and the solvent are replaced by spherically symmetric particles interacting with each other through the effective Lennard-Jones (LJ) potential. The utilization of the Weeks-Chandler-Andersen (WCA) procedure for the LJ potential and application of the fundamental measure theory [32] (FMT) to account for the short-range hard core correlations allowed the authors to fit the pressure and the surface tension of bulk water to the experimental values by changing the cutoff radius. The latter, in turn, enabled them to describe hydration of a set of hydrocarbons. Even though such cDFT approach remains a method of experimental data approximation, it has distinct advantages over the methods based on the implications of the PR EOS. The main one is that, in contrast to the PR EOS-based approach [33], which in its nature is a variant of the mean-field theory taking into account only long-range correlations of solution molecules, the cDFT-based method accounts for the short-range correlations of the solvent molecules due to the packing effects in the solute molecule solvation shell. Moreover, within the cDFT approach, all the interaction parameters are related directly to the effective LJ interaction potentials and, thereby, have a clear physical meaning.

The development, in a sense, of the cDFT application for the solvation free energy calculation can be traced by a number of papers [34], [35], [36], [37], proposing the so-called molecular density functional theory approach, which consists of several steps. The first one is computation of the site-site pair-distribution functions from the AAMD simulation of the pure solvent. The second one is obtaining of the direct correlation function from the Molecular Ornstein-Zernike integral equation, and, finally, calculation of the energy functional, the value of which, being the difference between the system grand thermodynamic potential with and without the solute molecule, gives the solvation free energy of the solute. Although such approach is faster than the full-scale classical AAMD simulation, the first step is still dependent on the force field choice. The results obtained for a series of halide ions are in good agreement with the experiment, but the question is what results could be achieved for drug-like molecules and how the calculation of such systems would influence the computational speed.

In this paper, we discuss two coarse-graining techniques that can be used for fast and sufficiently accurate estimation of solvation free energy values of sparingly soluble drug compounds in a scCO${}_2$ medium. One of the potential applications of the obtained data can be subsequent calculation of the solubility of the compounds, which is a prerequisite for using supercritical technologies, e.g. micronization or cocrystallization [38], [39], in order to enhance the final aqueous solubility and thus the bioavailability of drugs. The proposed approaches can make solvation free energy computation faster than in the aforementioned techniques, preserving the accuracy of the obtained values. The first one is based on cDFT, but, in contrast to the approaches discussed above, the only input values required to compute the interaction potential parameters are the critical parameters of the solvent and solute, which can be often found in literature. The second approach is based on the coarse-grained MD (CGMD) simulations of the LJ fluid, with the interaction parameters determined on the basis on the law of the corresponding states, knowing the critical parameters of the system compounds and critical parameters of the LJ fluid.