Challenges of the Use of Atomistic Simulations to Predict Solubilities of Drug-like Molecules

Background.Solubility is a physical property of extreme importance to the Pharmaceutical industry whose prediction for potential drugs has so far been a hard task.We attempted to predict the solubility of acetylsalicylic acid (ASA) by estimating 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 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 Einstein Molecule Method 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 experiment, 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. computational and it appears that both and approach applied the approach suitable for routine Abstract Background. Solubility is a physical property of extreme importance to the Pharmaceutical industry whose prediction for potential drugs has so far been a hard task. We attempted to predict the solubility of acetylsalicylic acid (ASA) by estimating 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 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 Einstein Molecule Method 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 experiment, 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 still is a challenge on the computational side, 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.


Introduction
Solubility is a critical property for pharmaceutical drug discovery, and problems with solubility can frustrate drug discovery efforts and block treatments. The bioavailability of a drug depends on the solubility difference between different crystal structures (polymorphs), dose, drug permeability, and formulation [1], so solubility plays a key role. Solubility problems can be unexpected and pose crucial obstacles that even threaten the administration of care. For example, a well-documented case occurred in the late 90's when ritonavir, an HIV-protease inhibitor marketed as Norvir, failed dissolution requirements [2]. Since ritonavir is not bioavailable in its solid form, it was administrated in capsules containing solutions designed not to be saturated with respect to the originally known molecular crystal (form I) [2]. The newly identified polymorph, form II, was unusually stable and unusually hard to crystallize; the preparation protocol of Norvir was inadequate to make capsules from the new polymorph, which severely threatened the supply of the drug in the market and endangered the lives of many HIV+ patients [2]. Considerable effort has already been devoted to the methods to predict crystal polymorphs [3,4,5,6,7,8,9], but much less attention has been given to methods to predict solubilities, with or without likely polymorphs as input.
Due to the importance of aqueous solubilities in different industrial processes and environmental applications, a scientific challenge consisting of the prediction of 32 solubilities given a database of 100 reliable measurements [10,11] was created with the goal of comparing the outcomes of different solubility prediction techniques. Participants employed methods such as artificial neural networks [12], quantitative structure-property relationship (QSPR) [13], and deep learning [14] to predict the aqueous solubilities of drug-like molecules. All of the employed methods were empirical and trained on existing measurements. The limitation of these methods, however, is the dependence of a training set of data that limits their applicability to compounds similar to those in the training of the set and impairs its transferability. Some newer methods attempt to predict solubilities based on a physical description of the interactions in solution and in the solid state, yielding results that are in principle rigorous given an accurate energy model and an adequate method. In these approaches, molecular systems are described by force fields, i.e, potential energy functions that contain parameters describing bonds, atoms, electrostatic and non-electrostatic interactions. Molecular dynamics (MD) or Monte Carlo (MC) simulations are commonly used to sample different configurations of the system described by an energy model called a force field. The simulations then allow the estimation of physical properties such as internal energy, free energy, and enthalpy under different conditions. The quality of the results of such methods depends on how well the force field describes the system under study and how good the sampling method is. Thus, some researchers have recently estimated aqueous solubilities using simulations of thermodynamic cycles encompassing the crystal, the ideal gas, and an infinitely dilute solution of a given molecule [15,16]. When the structure of the solid is unknown, some studies have substituted simulations of solid melts in place of a structure of the solid [17,18,19,20].
While these physical methods for predicting solubilities have received some attention in the literature, most are still in their infancy with only a handful studies applying them and it is not yet clear how broadly applicable they will be [17,18,19,20], and others have only been suggested or demonstrated in proof-of-principle tests [21,22,16,23]. Our view is that the time is ripe for physical methods to predict solubility, especially given the routine nature of solvation free energy calculations at present [24,25,26,27,28,29] which comprise essentially half of the solubility problem (see the Theory section). Polymorph and crystal structure prediction successes also mean that we may often have a suitable crystal structure of the compound as input [3,30,31,32,33,4,34,5,8,9,35], so what remains is to predict the solubility given a crystal structure and simulations of the relevant phases.
Here, we focus on adapting and testing an existing approach for solubility prediction in the hope that it will prove to be a generally applicable method for solubility prediction that can be applied routinely. This method uses all-atom molecular dynamics simulations to estimate absolute chemical potentials and predict aqueous solubilities of molecular solids, given the crystal structure (or an estimate thereof) as input.

Theory
The solubility of a molecular solid is related to the chemical potentials of each phase Solubility is defined as the maximum concentration of solute that can be dissolved in a selected bulk solvent. Chemical potentials (µ) of the solid-state solute and the solution are by definition equal at the solubility point, when the solution is in equilibrium with the solid.
Solid particles precipitate in concentrations higher than the solubility point because the solid phase becomes more stable in these conditions. In principle, we can predict at which concentration a molecule precipitates in solution if we calculate the chemical potentials of the components: where µ i is the chemical potential of component i; A is the Helmholtz free energy; G is the Gibbs free energy; N j, j =i is the number of molecules of each component in the mixture; V is the volume of the system; T its temperature; and P its pressure. Calculations from systems under constant V and T yield A; G is obtained from simulations under constant P and T conditions. In order to estimate the chemical potential of one component in solution and in its molecular solid, however, we need to know the absolute free energy of the system in these states. We calcu-lated absolute free energies using alchemical free energy calculations.

Alchemical free energy calculations can be used to calculate absolute free energies
The absolute free energy of a system can be determined if we know its partition function (Q), a function that connects microscopic properties of the system with macroscopic thermodynamic quantities. Unfortunately, it is very hard to calculate the absolute free energy of real systems because we don't know their partition functions. Free energy calculations allow us to bypass this problem, but require at least two states: a reference state whose free energy can be analytically or numerically found, and a final state of interest [36,37]. We chose to calculate the free energy difference using alchemical free energy calculations, a method in which we simulate a series of nonphysical intermediates between the end states [38].
Each intermediate state in the alchemical path is described by a Hamiltonian (q, p; λ), i.e, the energy of the state as a function of atomic positions (q), momenta (p) and a coupling parameter (λ): (3) where ini t ial and f inal respectively are the Hamiltonians of the initial and the final state; and f (λ) and g(λ) are functions used to mix the Hamiltonians and are usually set such that = ini t ial at λ = 0 and = f inal at λ = 1.
A variety of different estimators can be used to analyze alchemical free energy calculations, and have different strengths and weaknesses as well as different data requirements. Here, we employ several different estimators we introduce briefly in the following.
One way to calculate the free energy difference (∆A) between the end states is Thermodynamic Integration (TI) [39]: in which a set of discrete λ values correspond to states along the alchemical path. 〈〉 means that we are have to calculate the ensemble average of the derivative between the brackets. TI performs as well as more efficient methods if the integrand is smooth, but breaks down if this condition is not satisfied [40,41,42].
An alternate free energy estimation method computes ∆A directly via: where the ensemble average is calculated over the configurations of the initial state, and β is the reciprocal of k B T , the product between the Boltzmann constant and the absolute temperature. We call this approach exponential averaging [43] (EXP). Most free energy calculations involve many intermediates associated with the coupling parameter (λ), allowing simulation of intermediate states in between the two end states of interest. The free energy change between the end points of a path defined by N intermediates is: where ∆A n→n+1 is the free energy difference between (n + 1)-th and the n-th intermediate states. Eq. 5 can be used to calculate the free energy difference between each adjacent pair of states and yields the exact result at the limit of very large samples, but it is inefficient for most applications [38].
The Bennett acceptance ratio [44] (BAR) provides an estimator which is superior for most purposes. It calculates the free energy difference between the n-th and the (n + 1)-th states from the following relationship: N n e β(∆ n+1→n +∆A) n+1 (7) where N n and N n+1 respectively are the number of statistically independent samples in states n and n + 1, and ∆ n→n+1 = −∆ n+1→n are the Hamiltonian differences between n and n + 1. BAR is more efficient than EXP [45,46] and minimizes the free energy uncertainty [44]. Multistate Bennett acceptance ratio [41] (MBAR) is an extension of BAR that takes in consideration the degree of configuration space overlap between a given state and all other states in the transformation, while BAR only uses the information of neighboring states. MBAR and BAR perform similarly when the spacing between the intermediate states is moderate, but MBAR is the most wellperforming free energy estimator [42].

The absolute free energy of a solid is calculated using an ideal system as reference
In this work, we seek to predict solubilities of molecular solids. Part of this problem requires predicting the free energy or chemical potential of the solid. One way this has been attempted in the past is via the Einstein Crystal Method (ECM) which calculates the absolute free energy of a solid using an Einstein crystal as a reference state. In this method, the crystal lattice is made of atoms restrained to their positions by a harmonic potential; additionally, the center of mass of the system is held fixed [47].
In the ECM, and in this work, the absolute free energy of the molecular solid is found by designing a path where force field terms are progressively turned on, and the harmonic potential position restraints are turned off. The fixed center of mass is important to avoid a quasidivergence issue when calculating the free energy term of releasing the system from the harmonic position restraints, but the contribution of the fixed center of mass needs to be included in the cycle to obtain the correct absolute free energy for the system ( Fig. 1(a)) [47,48,49].
In ECM, the free energy is calculated by:

Figure 1. (a) Thermodynamic cycle representing the Einstein Crystal Method (ECM). (b) Thermodynamic cycle representing the Einstein Molecule Method (EMM). Notice that EMM requires only two free energy calculations despite being a bigger thermodynamic cycle. The canceling terms in (b) correspond to the free energies of fixing and releasing one atom in the crystal lattice [48].
where A EC F C M is the free energy of the Einstein crystal (EC) with a fixed center of mass (FCM); ∆A EC→I EC is the free energy difference between the Einstein crystal (EC) and the interacting Einstein crystal (IEC), i.e., the free energy difference in a transformation where the force field is progressively turned on throughout the calculation path. ∆A I EC→SF C M is the free energy difference between the IEC and the solid with a fixed center of mass (SFCM), i.e, turning off the harmonic restraints; and ∆A r elease C M is the free energy of release of the center of mass (CM).
ECM can be difficult to implement because of the need for a fixed center of mass, so our work here is based on an alternative approach which is easier to implement. When particles move in ECM, the lattice needs to be moved because the center of mass is fixed [47,50,48]. Our method of choice, the Einstein Molecule Method (EMM, see Fig. 1(b)), fixes a single atom in the lattice instead of the center of mass and is more easily implemented than ECM because of the relative difficulty of introducing center of mass restraints into existing simulation packages [50,22,48,51,52]. EMM has been used to predict phase diagrams of TIP4P and SPC/E water models [48], free energies of ice polymorphs, solid methanol and toy systems [49,52], and the solubilities of potassium and sodium chlorides [22,51].
In EMM, the free energy of a solid is: where A E M is the free energy of the ideal Einstein molecule; ∆A id→I E M is the free energy difference between the ideal Einstein molecule and the interacting Einstein molecule (i.e, turning on the force field); and ∆A I E M →solid is the free energy difference between the interacting Einstein molecule and the solid (i.e, turning off the harmonic restraints). The advantage of EMM over ECM is the ab-sence of the need to calculate a free energy term associated with releasing the fixed reference point [48].
Here, as per equation 9, we compute the free energy of the solid by combining the absolute free energy of the ideal Einstein molecule with two terms that we calculate via alchemical free energy calculations -∆A E M →I E M and ∆A I E M →solid ; these involve alchemically changing the interactions in the system. Numerical integration of Eq. 10 allows the calculation of the ideal term, A E M [52]: where A E M and Q E M are the free energy of the Einstein molecule and its partition function; U E M ,1 (Ω 1 ) is the potential energy of the fixed particle 1; U E M ,2 (r 2 , Ω 2 ) is the potential energy of a non-fixed particle at a distance r 2 of particle 1; Ω 1 and Ω 2 are all the possible orientations the molecules can have in the lattice; Λ, V , N , and β respectively are the de Broglie wavelength, the system's volume, its number of particles, and the reciprocal of k B T , the product of the Boltzmann constant and the absolute temperature.
The chemical potential of a component of a solution can be calculated using free energy calculations.
Another critical component of computing the solubility of a compound is estimating the chemical potential of a solute in solution, since the solubility point is the concentration at which the chemical potentials of compound in the two phases are equal.
The chemical potential of a component i in solution, µ i , has an ideal and and excess component: (11) where q i is the internal partition function of a single molecule of the solute, U(N i ) is the potential energy of the system with N i particles, Λ is the de Broglie thermal wavelength, and V is the system's volume [53]. 〈〉 ini t ial means that the term was obtained from an ensemble average over the configurations of the initial state (see Eq. 5). The first two terms of the equation above correspond to the ideal component of µ i ; the last one, µ ex i , corresponds to the excess component of µ i , and is associated with all non-ideal interactions of the extra component i with the solution (i.e. physical interactions that differ from those given by the ideal gas law). We obtained excess chemical potentials from solvation free energy calculations; the solute molecule is inserted in the solution by progressively turning on its interactions with the surrounding environment [24,54,28].
The challenge associated with the calculation of µ i is the calculation of the standard chemical potential of i, µ 0 i , the first term of equation 11. q i , the internal partition function, includes the rotation, vibrational, electronic and nuclear partition functions of a single molecule [53] and is unknown. Here, we found a way of calculating µ 0 i without the knowledge of q i by alchemically transforming a single solute molecule into a single Einstein molecule, whose absolute free energy we know how to calculate [50,48,49].

Distinctives of this work
We are aware of three main approaches to compute the solubility of solids in solution using physical approaches: ECM-based methods [21,23], EMM-based methods [22,51,55], and the approach of Michael Schnieders' and collaborators which computes sublimation and solvation free energies and uses these in an alternate thermodynamic cycle to obtain solubility estimates [15,56].
Many of the applications of these approaches have been to the solubility of ionic solids, with both ECM [21] and EMM-based approaches [22,51,55] having some success. However, molecular solids introduce substantial additional complexities for both of these approaches.
The ECM has seen an initial test on solubility estimation. Li et al. [23] used the ECM to estimate the solubility of napthalene, but made several approximations such as assuming that the internal partition function component of the solute cancels between environments (perhaps justified given napthalene's low solubility).
We are not aware of any work applying the EMM to solubility estimation of molecular solids; to our knowledge our work is the first to make such an attempt, though EMM has been used before to estimate the free energy of simple molecular solids [49,52] but not the solubility. This explains our need to find our own approach to estimate µ 0 i for a single solute molecule. The Schneiders approach is an orthogonal one that we do not examine here.

Systems under study
Here, we chose three systems to study: An argon crystal for some small initial tests, α-methanol to help establish our protocol, and acetylsalicylic acid (ASA) as our main object of study. ASA is a known anti-inflammatory whose most stable polymorph, form I [57], has an aqueous solubility of approximately 0.038 % mole fraction at 298 K [58]. We also used α-methanol at 150 K and a toy facecentered cubic (fcc) argon crystal [59] to help us find an optimal protocol to calculate the absolute free energy of a molecular solid. α-methanol was chosen because it had been used before in a study which applied the EMM to calculate the absolute free energy of the solid [52].
All simulations were run in GROMACS 4.6.7 [60,61,62,63]. With one exception, all simulations used the General Amber force field (GAFF) version 1.7 with AM1-BCC charges [64,65]; the exception was α-methanol, because we ran these simulations using the input files -coordinates and force field parameters -provided by Aragonès et al., who used an united atom version of the OPLS force field [52].
We simulated all solids and liquids using 5 ns Langevin dynamics simulations. ASA, α-methanol, and argon were simulated respectively at 298.15 K, 150.0 K, and 4.0 K. Our simulations had the same length as the simulations run by Aragonès et al. All solid state simulations were run in NVT conditions. Liquid state simulations were run in NPT conditions; pressure was kept constant at 101.335 kPa using the Parrinello-Rahman barostat [66]. We used the TIP3P water model [67] for all our liquid state simulations. More simulation details and example input files with full details can be found in the Supporting Information.

Calculation of the absolute free energy of molecular crystals
The absolute free energies of the solids were calculated from trajectories of simulation boxes with 64 ASA molecules, 100 OPLS methanol molecules, and 864 argon atoms with periodic boundary conditions. ASA's unit cell was obtained from Mercury CSD 3.8 [68] and the fcc argon crystal was obtained from the literature [59]. Simulation box sizes were chosen to be approximately between 2 nm and 3 nm to ensure that box sizes were large enough that atoms and their periodic copies were not within cutoff distance of one another. α-methanol's crystal was obtained from the Supporting information of Aragonès et al. [52] We used Amber14's ambertools [69,70,71,72] and ParmEd [73] to generate the ASA's and argon's solid state input files. All atoms but one were subjected to harmonic constraints in the x, y, and z coordinates. A single atom was kept fixed in space to act as the reference point for the calculations, as explained in the Introduction.
Monte Carlo integration yielded A E M , the free energy of the Einstein molecule, as it was previously done for α-methanol in the literature [52]. ∆A id→I E M and ∆A I E M →solid were estimated using Thermodynamic Integration (TI) [39] and the Multistate Bennett Accep-tance Ratio (MBAR) [74]. We used force constants of 4000 k B T/Å 2 to restrain atoms to their lattice positions in acetylsalicylic acid and argon simulations because it allowed us to use a reasonable time step of 1.0 fs in all simulations. α-methanol simulations used the same force constant that had been previously used by Aragonès et al. [52].
We used alchemical free energy calculations to obtain the difference in free energy between the reference Einstein molecule and the solid. This step was divided in two parts: (a) the force field parameters are alchemically turned on, and (b) the harmonic constraints are turned off.
Here, we deviate from earlier work which calculated the absolute free energy of a solid using EMM by introducing additional intermediate states to improve accuracy, along with using a superior free energy estimator.
For the calculation of ∆A id→I E M , we found it was crucial to introduce intermediate states; we also switched to using the MBAR estimator. The original EMM calculation of the absolute free energy of a solid [48,49,50,22,51,52] estimated ∆A id→I E M using exponential averaging (EXP) with just two states: the Einstein molecule (EM) and the interacting Einstein molecule (IEM) [48,49,50,22,51,52,21,55]. As EXP is known to have convergence issues and biases [40,38,41,45], we switched to the superior MBAR free energy estimator [74]. Additionally, when we did so, we found that overlap of states (as measured by the overlap matrix [75]) was insufficient so we created a series of intermediate states connecting both ends of the transformation.
For ∆A I E M →solid ., the original work used Thermodynamic Integration (TI) [39]. Here, we replaced TI with MBAR as our analysis method of choice. Generally, the literature shows that TI performs as well as more efficient methods like BAR and MBAR when the integrand is smooth [40,38,41], but it is sensitive to the choice and number of intermediate states [76]. MBAR is the most consistently well-performing free energy estimator [42] and exploits the overlap between states more thoroughly than its predecessor, the Bennett Acceptance Ratio (BAR) estimator [74]. Here, we chose to compare performance of MBAR and TI for calculation of ∆A I E M →solid for ASA and α-methanol; we also applied EXP as a comparison in the latter case only.

Chemical potential calculations
The chemical potential of a pure solid is its molar free energy: where N is the number of molecules in the solid, and A its Helmholtz free energy. The chemical potential of a substance i in water is defined as the derivative of the free energy of the system with respect to the composition: where G is the Gibbs free energy, and N i is the number of molecules of i in solution; P, T , and N H 2 O are the pressure, absolute temperature, and number of water molecules in solution, and are kept constant in the calculation. One important aspect to discuss is reason why we chose to calculate the Helmholtz free energy for the solid and Gibbs free energies for each solution. Solid state simulations with position restraints required running under constant temperature and constant volume conditions due to software limitations, therefore we were able to calculate A for the solids.. At constant pressure, both kinds of free energy are related by: Since solids are much less susceptible to volume changes than liquids, it is reasonable to consider that P∆V is negligible and ∆G ≈ ∆A. For instance, the difference in volume between the experimental ASA crystal structure and the simulation box after a constant pressure equilibration stage is 0.14 nm 3 . The P∆V term -i.e., the free energy difference discounting possible structure relaxation effects -would be much smaller than the simulation error.
As we explain in more detail in the Results section, successful absolute free energy calculations for molecular solids require a pathway involving a large number of alchemical intermediate states. The calculation of the absolute free energies of α-methanol at 150 K and acetylsalicylic acid required 600 states. Since GROMACS reads each λ until its fourth decimal place and the states need to be spaced more closely together as as the harmonic restraints are turned off (See Supporting Information), we decided to split each free energy calculation in sets of 100 states.
Liquid state simulation boxes were generated using the SolvationToolkit [77], a Python package that uses packmol [78], OpenMolTools (v0.6.7) [79] and OpenEye Python Toolkits [80,81,82]. Excess chemical potentials were obtained with the same solvation free energy protocol used in previous studies [28]: Starting from a fully interacting system, we progressively decouple the interactions of a single solute molecule with the remaining of the system, which allows us to calculate the free energy difference between a solute molecule in vacuum and in solution, i.e, the solvation free energy.
We also used alchemical free energy calculations using a single Einstein molecule as a reference state to estimate the standard chemical potential of a substance, µ 0 i : where µ F F o f f i and µ r est r aining i respectively are the chemical potential associated with turning off the force field and chemical potential of restraining the atoms of the molecule to their lattice positions. µ ideal i is calculated using the Monte Carlo integration procedure that we used to calculate A E M to a single molecule.

Chemical potential of molecular solids
The first step to predict aqueous solubilities with the aid of absolute free energy calculations was the assessment of the methodologies we chose to use. Since our method is the same one used by Aragonès et al. [52] and wanted to be sure that we could reproduce previous results, we ran simulations for α-methanol at 150 K and estimated the free energies of solids using MBAR. Turning off the harmonic restraints was the challenging step. Our MBAR is a free energy estimation method that minimizes the free energy variance and considers the overlap between a given state and all the others in the transformation path [41], which means that high uncertainties (±3k B T ) suggest the presence of problems in the transformation's path. TI's uncertainty estimates are much lower, but we believe this is an artifact. Error analysis for TI simply does not work the same way and does not give insight into whether exploration of phase space is adequate, unlike MBAR. Specifically, uncertainty estimates from TI usually factor in only the uncertainty in the integrand at each sampled lambda value and could potentially also factor in the smoothness of the integrand (i.e. numerical integration error) but do nothing to factor in whether the integrand will in fact vary smoothly in between lambda points; usually no data is available on this. BAR and MBAR, in contrast, factor in information about how well the intermediate states overlap in phase space and reflect high uncertainties when phase space overlap is poor. In our experience, usually TI would suffer from similar problems if additional intermediate states were added, but uncertainties in TI typically do not reflect this, as is the case here. Thus, the high uncertainty of the MBAR value indicates a sampling/convergence problem which warrants further exploration.
To explore the high uncertainty of our MBAR free energy estimates, we examined the degree of overlap the intermediate states had with each other. Phase space overlap analysis [83,84,85] quantifies the probability that any given configuration of an intermediate state can be found in other states. A good rule of thumb for designing a set of free energy calculations spanning between two states is to ensure that the states along the path have significant overlap with their neighbors as shown in Fig. 2. More overlap improves the quality of the MBAR free energy estimation: Fig. 2(b) represents a set of restraining simulations where the free energy uncertainty can potentially be accurately estimated using BAR and MBAR; Fig.  2(a) shows a case where it cannot. In our case we find that the α-methanol simulation using 18 intermediate states does not have adequate overlap (Fig. 3)-specifically, the states 4 ≤ λ i ≤ 17 do not have overlapping configurations with other states, which explains the 3 k B T uncertainty in our MBAR estimate.
Since prior work had appeared to do this estimation successfully [52], we were uncertain why we were encountering such overlap problems, so we studied an even simpler system. We calculated ∆A I E M →solid of fcc argon at 4 K with 18 states as in our α-methanol free energy estimation. MBAR yielded an error estimate of infinity while TI estimated ∆A I E M →solid to be −1666.5(8) k B T which, as we show below, is incorrect. This path resulted phase space overlap diagram without overlap between the states after state number 2 (Fig. 4). Apparently as the harmonic potential that holds atoms in their lattice positions tends to zero, atoms become rather mobile, dramatically decreasing phase space overlap and leading to poor free energy estimates.
To improve phase space overlap, we introduced more intermediate states along the path for removing the restraints (see Fig 2). We chose to break down the simulation in smaller parts, adding a significant amount of states near the point where the harmonic restraints are approximately zero. The MBAR estimate of ∆A I E M →solid for fcc argon is −1016.0(2) k B T using 300 states. TI's corresponding value was −1017(1) k B T, differing by far from the (incorrect) value of −1666.5(8) k B T obtained above with fewer states. Phase space overlap diagrams showed significant improvement in the configuration overlap between the states (Supporting Information). Thus, increasing the number of states was an effective strategy, and we used it in all subsequent calculations.
Even though our α-methanol results were similar to results previously published by other authors [52], we need to emphasize that reliable free energies resulted from simulations with a large number of intermediate states, as can be seen in Table 1. Despite its conceptual simplicity, calculating the components of the absolute free energy of a solid to a point where there is significant phase space . . . overlap between the intermediate states is computationally demanding. A 900-atom OPLS α-methanol system required 40 states to calculate ∆A id→I E M , and 600 states for ∆A I E M →solid . We chose these intermediate states in advance, and these ultimately led to free energy errors smaller than 0.1 k B T; the estimated TI and MBAR values differed by no more than 0.3 k B T. Our results for ASA using an optimal number of states can be seen in Table 2. The MBAR chemical potential of ASA at 298.15 K equals to −220.67(3) k B T. The computational cost of calculating A ASA was high; Each state required a separate simulation (of a 1344-atom ASA system), with 718 states in total. Simulations typi-cally required 11 hours on a single CPU, so the calculation of a single absolute free energy of a molecular solid required approximately 7898 CPU-hours.
Chemical potential of solutions and the solubility of GAFF acetylsalicylic acid in TIP3P water. Equation 11 states that the absolute chemical potential of a solution is determined by three quantities: µ 0 i , the standard chemical potential; µ ex i , the excess chemical potential of the component at a concentration of χ; and a volume-dependent ideal gas component of k B T · ln (Λ 3 i · N ASA /〈V 〉 solut ion ). µ 0 ASA only required information regarding the internal structure of the molecule [53], thus we estimated µ 0 ASA by alchemically transforming a single solute molecule into a single Einstein molecule (Table 3), whose absolute free energy we know how to calculate. We used the same number of states that we chose for the solid state simulations and we found that µ 0 ASA is equal to −150.7(2) k B T, as discussed in the last subsection of the Methods section.
Concentrations, volumes and excess chemical potentials can be seen in Table 4. We obtained the excess chemical potentials from solvation free energy calculations [24,54,28]. Volumes were obtained from the state in the alchemical path where the solute was fully coupled    to the rest of the system. The experimental aqueous solubility of ASA is approximately 0.038 % in water at 298 K [58], but our model predicts that ASA is effectively insoluble in water ( Figure 5). While all-atom simulations can yield solubility estimates given adequate simulation time and a correct method, the computed solubility will be that dictated by the underlying energy model or force field, and will not necessarily match experiment. Here, we use GAFF, a general-purpose force field with known limitations [86,71,87,28]; apparently, here, the right answer for the force field is not correct. Perhaps this is because of limitations in describing the solid state, as the force field is parameterized for liquid state simulations. Indeed, classical fixed charge force fields have shown severe limitations for polymorph prediction for these reasons [31,34,33,5,35] . Also, point partial atomic charges regularly used in molecular dynamics do not describe electrostatic interactions in a solid particularly well [88]. In the case of the ASA crystal, it is possible that its hydrogen bonds and π-stacking interactions add layers of complexity that are not properly described by GAFF.

Discussion and Conclusions
Despite its theoretical rigor, solubility prediction from absolute free energy calculations is a difficult task: it is computationally expensive and, at least in the present approach, requires many different steps and a great deal of care. Here, we attempted to develop and test a general approach to compute the solubility of molecular solids by adapting the Einstein Molecule Method (EMM) to tackle this problem, as discussed above.
To tune our methodology, we initially decided to repro- duce the absolute free energy of solid α-methanol, one of methanol's polymorphs, at 150 K using EMM before doing the same calculations for our compound of choice, ASA. We verified that the free energy differences between the Einstein molecule and the interactive Einstein molecule (∆A E M →I E M ) and between the latter state and the solid (∆A I E M →solid ) were more reliably estimated with the multistate Bennett acceptance ratio (MBAR). The absolute free energy of the crystal (as computed for united-atom OPLS α-methanol) agreed with results found in the literature, which suggested that we were on the right path. We did, however, require a very large number of intermediate alchemical states to obtain accurate free energy estimates, making these simulations fairly computationally demanding.
We then chose to calculate the solubility of ASA due to its pharmacological importance and due to its relative complexity in comparison to previous molecular solids whose absolute free energies have been computed via EMM previously [52]. As for α-methanol, this calculation required a large number of intermediate alchemical states and considerable computational cost -approximately 8000 CPU hours for a single absolute free energy calculation for the molecular solid, even with the crystal structure as input. Perhaps the number of intermediate states could be further optimized, but clearly a large number of intermediate simulations was required and thus considerable computational cost. Despite all of this, we still could not reproduce the experimental aqueous solubility of acetylsalicylic acid (ASA); experimentally it is modestly soluble, whereas our work would suggest it is essentially completely insoluble in water, likely due to force field limitations.
The solubility of naphthalene was recently estimated using a similar methodology, the Extended Einstein Crystal Method [23], but with additional approximations. Specifically, since naphthalene molecules interact very weakly with each other in the crystal lattice and with water molecules in solution, the differences between the internal partition function of a naphthalene molecule in the solid and in the solution were assumed to be negligible. This allowed the authors to drop some complexities in treatment of the solution-phase part of the calculation. However, that approach is only suitable for compounds that are only very weakly interacting in solution and in the crystal. ASA, in contrast, is a molecule that interacts strongly with other ASA molecules in its crystal lattice and with water molecules in solution via hydrogen bonds. For instance, an important crystalline feature that is not necessarily present in solution is the dimer structure, with two ASA molecules bound together via hydrogen bonds between the carboxylic acid groups. Differences between the internal partition functions of the molecule in the solid (q solid ASA ) and in solution (q solut ion ASA ) would probably not be negligible in this scenario, thus a more general approach is needed for treatment of such cases. Our work here provides one attempt in that direction.
Overall, the present approach seems to have significant limitations -most notably that the computational expense is considerable, and the resulting estimated solubility is quite inaccurate. Perhaps both of these may be surmountable; GPU-based free energy calculations can be dramatically faster potentially making an 8000 CPUhour calculation be 80 GPU hours which would amount to overnight on 8 GPUs, and perhaps this could be optimized via changes to simulation time and number of intermediate states. And with better force fields, perhaps accuracy could be improved; the AMOEBA-based approach of Schnieders shows considerable promise [15]. Alteratively, other approaches may be of interest. Solubility has been predicted by simulations using pseudocritical paths (i.e., paths were molecular crystals are transformed in tractable Einstein crystal-like states between the ending states of the transformation [89,90,91,92]) and a single experimental reference point [91]), and with the aid of a thermodynamic cycle formed by the molecular crystal, the molecule in vacuum, and the solvated molecule [15].
We believe the time has come for routine physical methods for estimation of solubility, even if improved force fields prove necessary before results have significant ac-curacy for application to biomolecular design problems.

Author contributions Competing interests
DLM is a member of the Scientific Advisory Board for OpenEye Scientific Software.

Grant information
DLM and GDRM appreciate the financial support from the National Science Foundation (CHE 1352608), and computing support from the UCI GreenPlanet cluster, supported in part by NSF Grant CHE-0840513. GDRM appreciates support from the Brazilian agency CAPES -Science without Borders program (BEX 3932-13-3)