Predicting Solubility of Small Molecules in Macromolecular Compounds for Nanomedicine Application from Atomistic Simulations

Solubility of small molecules in macromolecular compounds is of fundamental importance for a number of applications, including the growing field of nanomedicine. Here, approaches for in silico solubility predictions based on atomistic models are evaluated. A computationally efficient atomistic simulation procedure based on the concept of inherent structures is proposed for statistical sampling of polymer conformations. Comprehensive test simulations of several polymers and common solvents confirm the accuracy of the procedure for calculation of physicochemical properties such as cohesive energy densities, which along with the Hansen solubility parameters and the Flory–Huggins (FH) theory facilitate rapid, qualitative solubility predictions. However, the FH theory fails to model specific interactions such as hydrogen bonding. Therefore, more accurate predictions are obtained employing the perturbed hard sphere chain (PHSC) equation of state (EOS). As test case, aqueous poly(ethylene oxide) (PEO) solutions revealing strong hydrogen bonding are used. The physicochemical properties including the PEO‐water phase diagram calculated by the PHSC EOS show good agreement with experimental observations. Consequently, the combination of qualitative predictions using the FH theory for rapid prescreening with computationally more demanding parameterization of the PHSC EOS provides a promising tool for in silico solubility predictions with potential applications in nanomedicine.


Introduction
Predicting solubility of small molecules in macromolecular compounds is of fundamental importance for a number of DOI: 10.1002/adts.202000001 applications. One of them is the growing field of nanomedicine, [1][2][3] where achieving optimal solubility of small molecules in polymeric nanocarriers represents a prerequisite for an efficient targeted delivery of therapeutics or imaging agents. [4] However, experimental optimization of the encapsulation efficiency by adjusting the polymer structure and chemical composition using trial-and-error approaches is often costly and time-consuming. As a consequence, predictions of the thermodynamic compatibility between polymers and actives using computer simulations can speed up the process of discovery and optimization of new drug delivery systems. [5] In particular, atomistic simulations successfully predicted the thermodynamic drug-polymer compatibility and provide also a detailed understanding of the underlying intermolecular interactions (e.g., refs. [6][7][8][9]).
One of the most widely used quantity for characterization of the intermolecular interactions are the Hildebrand solubility parameters (SPs) defined as the square root of the cohesive energy density (CED) of pure chemical compounds. [10,11] In combination with the Flory-Huggins (FH) theory SP can be used to predict the enthalpy [9] and the free energy [12][13][14] of mixing as a measure for the thermodynamic compatibility between active substances and polymers. In addition, Hansen proposed an extension of the SP approach for more detailed analysis of the intermolecular interactions. [15] It is based on the separation of the CED into different energetic contributions such as dispersion and hydrogen bonding interactions. Recently, a machine learning procedure has been proposed for predictions of Hansen SP in order to provide a versatile tool for polymer solubility predictions. [16] Although several computational studies reported solubility predictions in agreement with experimental observations, [3][4][5][6] the general applicability and accuracy of atomistic simulations in combination with the FH theory has been debated. [17][18][19] For several mixtures of polymers and actives results obtained employing the FH theory parameterized by atomistic simulations were in disagreement with experimental observations. [18] This indicates the limitation of such solubility predictions connected with two key factors that determine reliability of thermodynamic modeling: i) accuracy of atomistic simulations and ii) derivation and parameterization of thermodynamic models.
Accurate atomistic simulations of amorphous bulk materials are in general a challenging task. In particular, the generation and sampling of statistically significant macromolecular conformations is indispensable for reliable modeling of amorphous polymers. Several studies used computationally very challenging long-timescale molecular dynamics (MD) simulations, with simulated times up to µs for relaxation and sampling of polymer conformations. [19,20] Configurationally biased Monte Carlo (MC) simulations have also been applied for generation of polymer conformations, which were subsequently equilibrated using MD simulations in the ns scale for sampling of the CED (e.g., refs. [6,7]). However, accuracy and reproducibility of such sampling approaches for the calculation of the CED have not yet been critically evaluated.
An alternative computational procedure initially proposed by Stillinger and Weber uses the concept of inherent structures (IS) for calculation of the atomic structure and thermodynamic quantities of amorphous bulk materials. [21][22][23] For a system containing N atoms, this approach separates the 3N-dimensional potential energy surface (PES) and the corresponding partition function into nonoverlapping basins. [22] Such basins are represented by IS, which are calculated by applying geometry optimizations to MD trajectories of the equilibrated liquid state above the glass transition temperature T g . This facilitates computationally efficient calculation of the configurational partition function. [24,25] Consequently, not only atomic structures but also thermodynamic quantities of glass forming liquids at temperatures close to T g can be calculated. [25] Simulations employing the concept of IS have been successfully used for elucidation of the structure and properties of various amorphous materials including polymers and biomolecules. [23,[26][27][28] However, to the best of our knowledge this methodology has not yet been applied to polymer solubility predictions.
The second key prerequisite for accurate predictions of the physicochemical properties of polymer mixtures with small molecules is the parameterization of thermodynamic models. While atomistic simulations explicitly consider intermolecular interactions including specific interactions such as hydrogen bonding, the FH theory fails to model such interactions due to the oversimplified mean-field lattice approximation. [17,18] This is due to the fact that the FH theory assumes positive entropy of mixing of the ideal lattice. However, strongly hydrogen-bonded solutions such as alcohol-water mixtures [29] are known to show negative entropy of mixing due to cluster formation in solution. [30] This applies also to aqueous solutions of poly(ethylene oxide) [31] (PEO) and PEO-drug mixtures [32] showing pronounced hydrogen bonding. Therefore, more sophisticated equation of state (EOS) theories are required for reliable predictions of the thermodynamic compatibility between the polymers and the active substances. [13] Among the numerous EOS theories for polymers [33] is the perturbed hard sphere chain (PHSC) EOS. [34][35][36] Parametrization of the PHSC EOS most frequently employed experimentally determined pressures as a function of the temperature as well as phase diagrams. Although few studies used atomistic simulations for derivation of the PHSC EOS parameters for polymers, [37,38] a critical evaluation of the accuracy of this approach for solubility predictions, in particular for mixtures showing pronounced specific interactions, is still lacking.
This work addresses the two essential factors necessary for reliable polymer solubility predictions-the accuracy of atomistic simulations for statistically significant sampling of polymer conformations and the thermodynamic modeling. First, a simulation procedure combining MC and MD simulations for sampling of IS of amorphous polymers is proposed. Comprehensive test simulations are performed in order to evaluate the accuracy and reproducibility of the statistical sampling of polymer conformations and related physicochemical properties such as CED and pressures. Subsequently, the Hansen SPs are calculated for qualitative solubility predictions of several polymer solvent combinations including novel polyesteramides (PEA) recently synthesized for nanomedical applications. [39] Finally, the accuracy of different thermodynamic models for calculation of the physicochemical properties is evaluated for aqueous PEO solutions as a test case for strongly hydrogen-bonded polymer mixtures. The models parameterized from atomistic simulations include the FH theory along with a composition-dependent correction of the interaction parameter as well as the PHSC EOS.

The Flory-Huggins Theory for Solubility Predictions
The FH theory [12][13][14] is based on the regular solution theory of nonelectrolytes [10] using a mean-field lattice approach that assumes randomly mixed and equally sized segments. In the following, the term segments is used for arbitrarily defined molecular fragments while the term unit refers to polymer repeating units and low molar mass molecules. In the FH theory, the average molar volume of the mixture v m is usually chosen as the size of one segment. Every segment pair interacts with an average energy ij , where i, j = 1 and 2 for polymer and active segments, respectively. The FH interaction parameter FH is given as [40] with z FH as the segment coordination number of the mean-field lattice. In our recent publication, [39] a composition-dependent FH interaction parameter z z = 1 RT Adv. Theory Simul. 2020, 3,2000001 www.advancedsciencenews.com www.advtheorysimul.com was derived from atomistic simulations for different PEA solutions. This concept denoted here as FH z considers the change of the intermolecular structure (using F 12 , F 1 (x 1 ), F 2 (x 2 )) with varying chemical composition implicitly including unequally sized segments (see the Supporting Information). [39] Similar composition-dependent FH parameters were derived from experiments to account for changing molecular environment of polymer chains in solution with varying solvent concentration. [41] A simplified expression for FH can be obtained by approximating CED of a mixture as C m = √ C 1 C 2 , where C 1 and C 2 are CED of the components. [11] Combined with the definition of the Hildebrand SP, i = √ C i , and SP of the mixture, 2 m = 1 2 , yields the Hildebrand-Scatchard equation for the energy of mixing (Equation (S2), Supporting Information) with the approximated FH interaction parameter HS where v s denotes the average molar volume of both components. In order to analyze the underlying interactions following the concept of Hansen, [15] the CED and SP can be separated into electrostatic (el) and van der Waals (vdw) contributions, i,el and i,vdw , respectively, with the total Hildebrand SP 2 i = 2 i,el + 2 i,vdw . [42] Applying the mixing rule to both components, 2 m,el = 1,el 2,el and 2 m,vdw = 1,vdw 2,vdw , yields This allows qualitative solubility predictions without explicit consideration of polymer-solvent mixtures (see the Supporting Information).

The Perturbed Hard Sphere (HS) Chain EOS
More accurate thermodynamic modeling compared to lattice fluid models such as the FH theory can be achieved using polymer EOS. [33,43] A pressure explicit EOS describes the dependence of the pressure p on the temperature T and volume V or number density , respectively, and is frequently expressed using the dimensionless compressibility factor Z [43] In case of the ideal gas, the compressibility factor is simply Z id = 1. In this work, the PHSC EOS [34][35][36] is parameterized based on atomistic simulations. The model contains corrections to Z id considering pairwise contacts of HS segments Z HS , chain bonding (CB) of spherical segments Z CB as well as attractive intermolecular interactions by using the perturbation (pert) term of the van der Waals EOS Z pert More detailed description of the PHSC EOS including the exact expressions for Z is given in the Supporting Information and in refs. [34][35][36]. The PHSC EOS requires three adjustable parameters for every segment pair i − j with i, j = 1 and 2. Namely, the segment diameter (center-to-center distances) ij , the number of HS segments per molecule r i , and the (attractive) intermolecular interaction energy ij . These parameters are most frequently derived for each component by fitting Equation (6) to observed p(T, ) curves and for the parameters between unlike segments by using experimentally determined phase diagrams. [44] However, the PHSC EOS fails to model polymer mixtures revealing pronounced specific (directional) interactions such as hydrogen bonding between the HS segments. [35] As a consequence, lower critical solutions temperatures (LCSTs) or closedloop phase diagrams as in case of aqueous PEO solutions cannot be determined. One method to include the effects of specific interactions in EOS and lattice fluid theories such as FH is to define temperature-dependent interaction parameters ij (T) and ij (T). [45] A simple model for ij (T) involves two potentials for the intersegment interaction, one considering specific (directional) s,ij and one for nonspecific ns,ij interactions. [34,[45][46][47] The interaction energy of segment pairs at temperature T is given by the average interaction parameter with ij = s,ij − ns,ij > 0 and the fraction of specific interactions f ij (T). This fraction is calculated as the ratio between the segments interacting by nonspecific and specific interactions assuming the Boltzmann distribution, yielding [45,47] with the degeneracy ratio q ij of nonspecific to specific interactions. In combination with lattice fluid models, this approach has been successfully applied to model closed-loop phase diagrams of PEO-H 2 O mixtures. [47] In this work, this model is applied not only to derive ij (T) but also to calculate the FH parameters FH and z (Equations (1) and (2)). Parameterization of the PHSC EOS provides the expression for pressure as the function of temperature and number density p(T, ), which allows calculation of all thermodynamic quantities. The excess internal energy calculated using pressure explicit EOS of a liquid with number density 0 , containing mole fractions x i of molecules i is given as [43] e EOS ( which is related to the CED, C = − 0 e EOS . Similarly, the excess Helmholtz free energy can be calculated as [43] a EOS ( In this work, e EOS and a EOS were calculated numerically for selected T, 0 , and x i . In addition, the enthalpy h EOS = e EOS + pV Adv. Theory Simul. 2020, 3, 2000001 Figure 1. Computational workflow used for parameterization of the FH theory (FH and FH z ) and the PHSC equation of state employing the SW and LJ intersegment potentials ij (R), respectively. Amorphous structure models generated by MC and MD simulations are used for calculation of the RDF of the coarse-grained models as well as cohesive energy densities C(T) and pressures p(T) as a function of temperature. and the Gibbs energies g EOS = a EOS + pV were determined. In order to test the accuracy of the PHSC EOS derived from simulations, the water vapor pressures as a function of temperature as well as the PEO-H 2 O phase diagrams were calculated and compared to experimental data. However, the volume of incompressible phases depends only slightly on pressure and thus the volume work of mixing pdV at moderate pressures considered in this work (below 200 bar) is negligibly small. Therefore, the vapor pressure of pure water was taken for calculation of the PEO-H 2 O phase diagrams at all mixture compositions. All phase equilibrium calculations used the conditions of the mechanical p ′ (T, ) = p ′′ (T, ) and the chemical equilibrium ′ i = ′′ i of component i for the phases ′ and ′′ .

Model Parameterization by Atomistic Simulations
The required parameters of the thermodynamic models described in Sections 2.1 and 2.2 were derived from atomistic simulations using the computational workflow depicted in Figure 1. A configuration-biased MC algorithm [48,49] was used to construct initial amorphous structure models. Next, the CED C(T) and pressures p(T) as a function of temperature were derived from MD simulations. Calculations of the intermolecular radial distribution function (RDF) g ij (R) as a function of distance R used coarse-grained models generated with the Materials Studio program suite [50] for definition of the molecular units employing atomistic structure models obtained from MD simulations. The CED is related to the RDF through the intersegment potential ij (R) [51,52] for a system containing l different components.
In case of the FH theory, ij (R) was approximated as squarewell (SW) potential with the temperature-dependent interaction energy ij (see Equation (1)). Using the definition of the segment coordination number z ij calculated employing the RDF of the coarse-grained models, Equation (11) reduces to a simpler form, e.g., C = −0.5 z ij ij for a pure substance. [39] Calculations of z ij (Table S4, Supporting Information) used the cutoff distance R 1 that corresponds to the diameter of a sphere with the volume equal to the average volume of one segment. In this way, the SW pair potential can be used for derivation of the parameter ij required for the FH theory. Further details on the calculation of the FH interaction parameter and derivation of its compositiondependent correction (FH z ) for consideration of unequally sized segments can be found in our recent publication (ref. [39]).
The parameterization of the PHSC EOS used the Lennard-Jones (LJ) potential for the interaction between the molecular units (un) along with the parameters un ij and un ij (Table S5, Supporting Information). For un ij (T), the temperature-dependent model given by Equation (7) is applied along with the adjustable parameters un ns,ij , un ij , and q un ij . In addition, the parameter un ij is used that corresponds to the equilibrium (center-to-center) distance of two units. The parameters were derived based on fitting of the CED (using Equation (11)) and the pressure p as a function of temperature T, both derived from atomistic simulations. For this purpose, the compressibility factor (Equation (5)) defined for the intermolecular LJ potential was used [51] where ′ ij is the first derivative of the effective LJ potential with respect to the distance R. The cutoff distance R 1 in Equations (11) and (12) was set to 15 Å. Next, for derivation of the PHSC EOS parameters ( ij , r i , ns,ij , ij , and q ij ), the LJ potential ij (R) derived from atomistic simulations is separated using the Weeks-Chandler-Andersen decomposition. [53] This decomposition was also applied for derivation of the PHSC EOS, [35] in which ij is split into a repulsive ij,rep and an attractive ij,att part, such that ij = ij,rep + ij,att . The reference part of the PHSC EOS Z ref = 1 + Z HS + Z CB depends only on the repulsive potential ij,rep , www.advancedsciencenews.com www.advtheorysimul.com while Z pert contains only attractive interactions ij,att . [35] Consequently, the pressure can be split into the repulsive and attractive parts p att (using Z pert ) and p rep (using Z ref ), respectively. The excess potential energy e EOS (Equation (9)) and, thus, the CED can be split in the same way yielding C att and C rep . This allows for targeted parameterization of the reference part Z ref of the PHSC EOS for noninteracting HS chains and the perturbation part Z pert accounting for all attractive interactions. The same decomposition of the CED and pressures was applied here employing the effective LJ potential obtained from simulations along with Equations (11) and (12). It was used for the derivation of the PHSC EOS parameters for PEO and H 2 O by least-square fitting of p att and p rep as well as C att and C rep in a temperature range from 300 to 550 K (300 to 400 K in case of polylactic acid (PLA)). In addition, the number density at 300 K was included in the least-square fitting procedure. Since the CED for PEO was considerably overestimated by the PHSC EOS, only the pressures p att and p rep were used for fitting of the parameters of unlike segments (see the Supporting Information) in case of the PEO-H 2 O mixture.

Details of Atomistic Simulations
All atomistic simulations were performed employing modules of the Materials Studio (Version 17.1) [50] program suite along with the COMPASSII force field. [42] MD simulations employed the Forcite program module along with the Nosé-Hoover thermostat [54,55] and a time step of 1 fs. Construction of amorphous, 3D unit cells used the configuration-biased MC method implemented in the Amorphous Cell module, which is based on the Theodorou and Suter [48,49] algorithm and the Meirovitch scanning method. [56] The MC generated structure models were geometrically optimized, followed by a simulated annealing procedure. For this purpose, systems were equilibrated at T = 300 K using MD simulations along with the canonical ensemble (NVT ensemble). Subsequently, the temperature was stepwise increased up to 1000 K and subsequently decreased back to 300 K. In each step, the temperature was increased (decreased) by 100 K followed by an equilibration for 5 ps. Next, the refined structure models were equilibrated for 100 ps using the isothermalisobaric ensemble (NPT ensemble) at zero target pressure and T = 300 K employing the Berendsen barostat [57] and further simulations using the Parrinello-Rahman barostat [58] with a duration of 300 ps. Average cell parameters were evaluated for the last 200 ps of the NPT simulations. Finally, the unit cells were scaled to the average cell parameters obtained from NPT simulations and the structure models were again equilibrated for 250 ps employing the NVT ensemble with a target temperature of 300 K. Average properties such as CED were calculated from the last 200 ps of the NVT simulations.
Sampling of the polymer conformations used two different simulation procedures denoted here as direct sampling (DS) and IS sampling. The DS uses one single MD simulation of an MC constructed structure model containing a sufficiently large number of polymer molecules. In order to determine the optimum system size for DS, this procedure was applied to PEO with a degree of polymerization (DP) of 25 and 50 as well as PLA with a DP of 25 using structure models containing 1-30 molecules in a simulation cell. To test the reproducibility of DS for calculation of the CED, the whole computational procedure starting from the MC construction of the initial models was repeated ten times for every system size. In addition, to evaluate accuracy of the COM-PASSII force field the DS was applied to several commonly used solvents listed in Table S1 in the Supporting Information using simulation cells containing 800 molecules.
In case of PEO, PLA, H 2 O as well as the PEO-H 2 O and PLA-H 2 O mixtures, the IS sampling was applied. A simplified representation of this approach is depicted in Figure S2 in the Supporting Information. First, atomistic structure models of the amorphous bulk polymer showing distinct polymer conformations (MC1-MC3 in Figure S2, Supporting Information) were generated using the configuration-biased MC algorithm. Subsequently, MD simulations were performed for every MC constructed structure model as described above. Calculation of pressure-temperature curves used additional MD simulations for temperatures ranging from 350 to 550 in 50 K steps with a duration of 250 ps employing the NVT ensemble along with unit cell volumes obtained from the NPT simulations at 300 K. The sampling of IS was performed applying structure optimizations every 1 ps to the MD trajectories calculated for 400 K yielding their potential energy per atom e IS . This temperature is slightly above the glass transition temperature T g of PEO (250 K) [59] and PLA (330 K). [60] Since the MD simulations used short simulation times and moderate temperatures close to T g , the IS sampling is restricted to the basins of the PES corresponding to particular polymer conformations (e.g., MC1-MC3 in Figure S2, Supporting Information). That is, the kinetic energy is too low to overcome higher energy barriers of the PES yielding broken ergodicity. [61] Even above T g , MD simulation with very long equilibration times up to µs are required due to the relatively slow relaxation of polymer backbone conformations. [19,20] However, assuming local (internal) ergodicity within the PES basins [62] the IS sampling yields the distribution P i (e IS ) for a certain polymer conformation i (e.g., MC1 in Figure S2, Supporting Information). The total (average) probability distribution P(e IS ) that an IS with energy e IS characterizes the macroscopic liquid state was calculated for n MC polymer conformations sampled by the IS as using the Boltzmann weights w i = exp(− e IS,i )Q −1 along with Q = ∑ n MC i exp(− e IS,i ) and the average IS energy e IS,i of a conformation i. Average quantitiesȲ, such as CED or pressure required for model parameterization (see Section 2.3.1), were calculated using the w i and the average properties calculated from MD simulations Y i according toȲ = ∑ n MC i w i Y i . The average intermolecular RDF of molecular units g(R) = ∑ N i P i g i (R) was calculated using the RDFs of the coarse-grained models g i (R) of N IS. One coarse grained bead was used to approximate each polymer repeating unit or a solvent molecule. Construction of coarsegrained models used the module Mesostructure Builder of the Materials Studio program suite. [50] The major advantage of the IS sampling is a trivial parallelization over n MC independent MD simulations. Several test simulations for IS sampling were performed, in particular for

Direct Sampling
A sufficiently large set of points in phase space (atomic configurations and momenta) has to be sampled by atomistic simulations to yield reliable values of the macroscopic, thermodynamic quantities. In order to find the best compromise between computational effort and accuracy, the reproducibility of the evaluation of CED expressed as the Hildebrand SP was tested by repeating the DS for ten different system sizes, i.e., number of polymer repeating units per simulation cell. The resulting SP for PEO (DP 25 and 50) and PLA (DP 25) as a function of the system size are shown in Figure 2a.
The calculated SPs showed large fluctuations (standard deviations) for PEO and PLA models containing less than 150 polymer repeating units per simulation cell. In contrast, models with 500 repeating units yielded standard deviation for SPs of less than 0.  [63] and PLA, [64] respectively. This demonstrates the reliability and robustness of the simulation procedure and the COM-PASSII force field. [42] To further verify the accuracy of the force field, the DS was applied to predict SP of various commonly used solvents. Figure 2b shows the comparison of the calculated and experimental [15] SP (see Table S1, Supporting Information). For compounds containing only C, N, O, and H, the calculated SP revealed a very good agreement with the experimental data. The mean absolute error (MAE) for this set of solvents was 0.6 MPa 0.5 and 30 J cm −3 for the SP and CED, respectively. However, considerable deviations of 4.2, 3.7, and 8.7 MPa 0.5 from the experimental values were obtained for the halogen-containing solvents hexafluoro-2-propanol (HfiP) and dichloromethane (DCM) as well as the sulfur compound dimethyl sulfoxide (DMSO), respectively. The COMPAS-SII force field was chosen since its parameterization included experimentally observed vaporization enthalpies. [42] To further verify the force field accuracy, the SP of ten molecular crystals were calculated including cyclic and noncyclic (C, N, H, O containing) compounds. Table S3 in the Supporting Infor- mation summarizes the SP calculated using the COMPASSII, Dreiding, [65] and PCFF [66] force field as well as the SP calculated using experimentally observed sublimation enthalpies [67] (Equation (S1), Supporting Information). COMPASSII shows the lowest MAE of 1.3 MPa 0.5 while Dreiding and PCFF demonstrate an MAE of 4.2 and 2.1 MPa 0.5 , respectively. We conclude that quantitative solubility predictions using the COMPASSII force field are only possible for polymers and active substances consisting of C, N, H, O, which is sufficient for many pharmaceutical and biomedical applications.

Inherent Structure Sampling
Applying DS to complex, amorphous polymer-solvent mixtures can be very challenging since large simulation cells are required to ensure the reliable sampling of all relevant molecular configurations. Therefore, IS sampling was investigated as potential, computationally more efficient alternative to DS. It employed a set of considerably smaller simulation cells containing 150 polymer repeating units. Figure S3 in the Supporting Information depicts the total probability densities P(e IS ) of the IS calculated for pure PEO and PLA as a function of their potential energy (per atom) e IS .

www.advancedsciencenews.com www.advtheorysimul.com
Here, e IS of an IS was calculated as energy difference between each IS and the lowest energy structure obtained during the IS sampling. In addition, the weighted distributions w i P i (e IS ) for each simulation cell are shown. In case of PEO, similar P i (e IS ) were obtained independent of the initial structure, due to a high flexibility of PEO backbone. In contrast, for PLA the mean value and sampled range of e IS showed strong dependence on the initial configuration, despite a sampling temperature of 400 K that is clearly above the T g of PLA (330 K). [60] Therefore, not all relevant polymer conformations could be sampled during an MD simulation of one single MC-generated structure using simulation times of less than 1 ns and such small unit cells (150 units). However, the rather short MD simulations are used to sample only the local atomic vibrations of the polymer configurations with relaxation times of few tens of ps. An improved sampling of the configuration space for a single polymer structure model would require MD simulation times in the order of µs. Yet, the IS sampling employs several (n MC = 10-20) MC constructed polymer configurations to obtain a statistically relevant sampling. In case of pure water, almost identical P i (e IS ) were calculated (not shown) due to the far shorter relaxation times of low molar mass compounds in the liquid state compared to macromolecular compounds. Figure S4 in the Supporting Information shows the Hildebrand SP and segment coordination number, z 11 , both derived from the intermolecular RDF as well as the interaction parameter 11 of the SW potential (see Figure 1) calculated for PEO and PLA using DS and IS sampling with n MC simulation cells. Starting from n MC = 10, the SP of PEO and PLA showed very good agreement with the values obtained from the DS. Moreover, the variation of SP for different n MC is significantly smaller than the intrinsic error of 0.6 MPa 0.5 estimated for the COMPASSII force field (shaded area in Figure S4a, Supporting Information). Similarly, the coordination number z 11 quantifying the intermolecular structure of the first coordination shell of the polymer repeating units agreed well with the results of the DS, albeit the IS sampling slightly underestimated z 11 by about 0.1 in case of PLA. The interaction parameter 11 revealed the same trend. It can be concluded that IS sampling using n MC = 10 provides statistically relevant sampling of both, the intermolecular structure and interactions, with same accuracy as the DS.
A more detailed comparison between IS sampling and DS is displayed in Figure S5 in the Supporting Information for the RDF of PEO repeating units and H 2 O molecules in the pure and mixed states. A very good agreement was obtained in each case. In fact, IS sampling provided higher resolution of the RDF than DS indicating a more accurate sampling of the intermolecular structure. Figure 3a shows the plot of electrostatic and van der Waals contributions to Hildebrand SP, el and vdw , respectively, for PEO and selected solvents (see Tables S1 and S2, Supporting Information), with each compound assigned to a point ( i,vdw , i,el ). Since the distance between two points is proportional to the energy of mixing Δe mix,FH (see Equations (S2) and (S4), Supporting Information), one can define a critical distance R c (Equation (S6), Sup-  Table S1, Supporting Information) and a PEO chain length of 50, R c is 4.6 MPa 0.5 . As a consequence, a circle can be defined around PEO ( el = 18.5 MPa 0.5 , vdw = 7.9 MPa 0.5 ) enclosing all solvents that are predicted to solubilize PEO for every solution composition. Figure 3a demonstrates that these qualitative solubility predictions yielded good agreement with experimental data for solvents with el < 15 MPa 0.5 . However, they failed for solvents with higher el , in particular for water ( el = 46.6 MPa 0.5 ) and methanol ( el = 26.9 MPa 0.5 ). This is connected with the formation of hydrogen bonds in solution. However, for sufficiently low el , SPs allow rapid qualitative solubility predictions. In addition, el can be used as an indicator for potentially important specific interactions, which require more advanced simulation techniques in order to achieve accurate predictions of the thermodynamic compatibility. Figure 3b shows the relative SPs Δ el and Δ vdw for various polymer-solvent combinations (see Tables S1 and S2, Supporting Information). This plot is divided into areas I and II by the line at Δ el = 10 MPa 0.5 . For Δ el > 10 MPa 0.5 , all polymers are insoluble in the solvents tested, except for PEO-water, PEO-ethanol, and PEO-methanol. That is, if el of a solvent (or drug) exceeds el of the polymer more than 10 MPa 0.5 , then they are expected to be immiscible, unless they mix due to strong specific interactions. Figure 3c displays the enlarged lower part of Figure 3b with II further separated such that all solvents that solubilize the respective polymer are located within area IIb with |Δ vdw | < 3 MPa 0.5 , while in IIa contains only immiscible polymer-solvents combinations. However, also insoluble polymer-solvent systems are present in IIb. This indicates that not only specific interactions play a vital role for accurate solubility predictions but also steric effects of polymer mixtures, which can be modeled using the compositiondependent FH parameter employing simulations of the actual mixture. [39]

Refined Modeling of Polymer Mixtures
In case of mixtures with strong specific interactions, e.g., aqueous PEO solutions involving strong hydrogen bonding, the simple solubility model described in the previous section clearly fails. The specific interactions can only be accounted for by atomistic simulations of the actual mixtures. However, even if the simulations of mixtures are performed, the simple thermodynamic models such as FH mean-field lattice theory do not consider the effects of specific (directional) interactions. Similarly, the PHSC EOS assumes a statistical ensemble of randomly mixed HS chains [34][35][36] without explicit consideration of specific interactions. In this study, the effects of specific interactions on the temperature dependence of the interaction parameters ij and ij of the FH theory and the PHSC EOS, respectively, were modeled using Equation (7) for calculating the closed-loop phase diagram of PEO-H 2 O. Figure 4 shows (total) CED as well as pressure p(T) of PEO, H 2 O, and PEO-H 2 O mixture as a function of temperature calculated using the parameter sets P1-P3 of the PHSC EOS (see Tables S6 and S7, Supporting Information) and the LJ potential (see Equations (11) and (12)) derived from atomistic simulations (SIM). Calculations used constant (zero pressure) density at 300 K obtained from NPT simulations and P1-P3, respectively. In case of pure PEO, the p(T) curves revealed a fairly good agreement at 300 and 350 K. However, due to the clearly lower slope of p(T) for P1-P3 compared to SIM, the PHSC EOS considerably underestimated the pressures at higher temperatures. The corresponding CED calculated using P1-P3 was systematically higher than SIM by about 30%. In contrast, for H 2 O, a larger slope of p(T) was found for P1-P3 in comparison with SIM, while the CED was up to 15% lower in case of P2 (at 300 K). Due to the systematic overestimation of the CED of PEO, the CED of the mixture was not included in the training data for fitting of the PEO-H 2 O parameters. Therefore, the p(T) curves calculated using P1-P3 demonstrated a very good agreement with SIM, and similar to PEO, the CED was overestimated up to 27%. The separation of the intersegment (PHSC EOS) and LJ potential into repulsive (rep) and attractive (att) parts allows a split of the total CED and pressure into attractive and repulsive contributions C rep , C att and p rep , p att , respectively. Deviations of C rep , C att , p rep , and p att calculated using the fitted PHSC parameter set P1 with respect to the results of the LJ potential derived from SIM at different temperatures are summarized in Table 1. For pure PEO, only minor deviations were obtained for C rep , while C att was considerably overestimated by P1 over the whole temperature range yielding too high total CED (Figure 4). In contrast, the too low values of p rep were compensated by too large p att at lower temperatures. However, deviations of p rep considerably increased with temperature, while the error of p att remained almost constant resulting in a lower slope of p(T) for P1 compared to SIM. Similarly, the relatively small deviations of p rep and p att in case of H 2 O compensated at lower T, yet the overestimation of p rep by P1 strongly increased at higher temperatures. Conversely, deviations C rep were virtually independent of temperature, while the overestimation of C att increased with increasing temperature yielding error compensation at higher T. In case of the PEO-H 2 O mixture, the largest deviations were obtained for C att and comparatively pronounced relative errors for C rep , both almost independent of T.
In summary, an accurate fit of the PHSC EOS to the results obtained from atomistic simulations of both, p(T) dependence and absolute values of CED for PEO containing systems could not be achieved. As shown by Equation (9), the (residual) potential energy e EOS < 0 calculated from a pressure explicit EOS decreases with increasing slope of p(T) at constant density. The steeper increase of pressure with temperature results in higher CED since C = − e EOS . Therefore, the PHSC EOS in this form is not accurate enough to reproduce the simulation results for temperature dependence of the CED and pressure with equal accuracy. This is mainly connected with two major shortcomings of the parameterized EOS. First, the simple van der Waals perturbation Z pert (Equation (6) and Equation (S12), Supporting Information) does not consider hard sphere contacts of a dense liquid, i.e., does not include a realistic (hard sphere) RDF [68,69] g ij ( , ij ) (Equation (S8), Supporting Information) as used in Z HS . More accurate expressions for Z pert have been derived by Hino and Prausnitz [70] using an SW potential for attractive intersegment interactions along with an approximate expression for consideration of g ij ( , ij ). Second deficiency of PHSC EOS is the simplified model for consideration of specific, directional interactions using an average, temperature-dependent interaction parameter ij (T) (Equation (7)). Obviously, this model is insufficient to accurately include the effects of hydrogen bonding. This can be resolved by adding analytical expressions for association sites of hydrogen bonds to the PHSC using, e.g., the Statistical Association Fluid Theory (SAFT). [71][72][73] Employing such extensions of the PHSC EOS appears promising to facilitate more accurate modeling of the temperature dependence of CED and pressures obtained from atomistic simulations. Figure 5 shows experimentally determined water vapor pressures [74] and phase diagrams [75] for PEO-H 2 O mixtures in comparison with the results calculated using the PHSC EOS and parameter sets P1-P3. The set P2 revealed the best agreement with experimentally determined water vapor pressures, although it deviated most from CED obtained from simulations ( Figure 4). However, also P1 and P3 provided good agreement with experiments up to temperatures of about 500 K. The calculated phase diagrams using the sets P1 and P2 depicted in Figure 5c show fairly good agreement with the experimental one. In particular, P1 accurately reproduced the experimentally observed upper critical solution temperatures (UCST) up to PEO weight fractions of 0.3, yet clearly overestimated the LCSTs. In case of P3, PEO was predicted to be soluble in water over the whole temperature (300-700 K) and weight fraction (0-1) range. In addition, the calculated phase diagrams using P1 for different PEO chain lengths showed a very good qualitative correlation with the chain length dependencies observed in experiments, particularly for the UCST. This clearly shows that our parameterization procedure of the PHSC EOS by atomistic simulations is capable of providing reasonable thermodynamic models even for mixtures showing pronounced specific interactions. Although the parameter sets P1-P3 are very similar providing almost identical fits of CED and pressures (with exceptions of the water CED for P2), they yield distinct results for the physicochemical properties of PEO-H 2 O mixtures. In particular, only small deviations of the Gibbs energies of mixing Δg mix (per molecular unit) can lead to very different phase diagrams. Figure 6a depicts Δg mix,FH and Δg mix,z calculated using the FH theory with a constant FH (Equation (1)) and compositiondependent z (FH z , Equation (2)) interaction parameter, respectively. In addition, values for Δg mix,EOS calculated using the PHSC EOS with parameter sets P1-P3 are shown. The parameters of the SW potential and the coordination numbers derived from simulations (see Table S4, Supporting Information) yielded positive FH and, therefore, incorrectly predicted positive values for Δg mix,FH for all mixture compositions. One shortcoming of this approach is the assumption of a random mean-field lattice with a coordination number z FH independent of composition, which is in clear contrast to the RDF calculated from simulations (see Figure S5, Supporting Information). The PEO-H 2 O RDF and the structure model depicted in Figure S5f in the Supporting Information show that PEO and H 2 O are not randomly mixed but rather partially separated at the molecular level. Such incomplete mixing at the molecular scale was observed for strongly hydrogen-bonded mixtures such as concentrated alcohol-H 2 O mixtures. [29] As second test case, the PHSC EOS was parameterized for PLA-H 2 O (see Tables S6-S8, Supporting Information). In agreement with the experimentally observed immiscibility of PLA and water, the calculated Δg mix are positive for virtually all compositions (see Figure S6, Supporting Information). Only at very low water concentrations of about 0.3 wt% the mixture is slightly exergonic with Δg mix = −50 J mol −1 (at 300 K). This coincides with water absorption experiments demonstrating a water uptake of PLA of about 1 wt%. [76] The simple polynomial fit of the coordination numbers z ij (see ref. [39]) for calculation of z provided a qualitative correction of Δg mix,FH for PEO-water mixtures below the PEO weight fraction  Figure 6b shows the enthalpic Δh mix and entropic contributions −TΔs mix to Δg mix for PHSC EOS with the parameter set P1 and FH z along with experimentally determined [77] Δh mix for a PEO chain length of 68 at 353 K. The composition dependence of Δh mix,z is similar to Δg mix,z due to the minor (negative) contributions of −TΔs mix,z using the ideal entropy of mixing Δs mix,z > 0. In contrast, positive values −TΔs mix,EOS up to a PEO weight fraction of about 0.8 were obtained in case of P1, while Δh mix,EOS , contrary to Δh mix,z , revealed fairly good agreement with the experimentally determined values. As for Δg mix,EOS , similar values of the enthalpy and entropy of mixing were calculated for P2 and P3. Table 2 compares Δg mix , Δh mix , and Δs mix calculated using the PHSC EOS with the parameter sets P1-P3 as well as the FH theory (FH and FH z ) with values for Δh mix obtained from simulations and experiments [77] for a PEO weight fraction of about 0.8 and T = 353 K. The error of Δh mix,sim ≈ Δe mix,sim calculated using atomistic simulations, was estimated using error propagation of the general definition of Δe mix = e m − x 1 e 1 − x 2 e 2 assuming an uncertainty of the CED for pure PEO, water as well as their mixture of 30 J cm −3 yielding error limits of about ± 1.9 kJ mol −1 . Therefore, the incorrect simulation result of an endothermic mixture with Δh mix,sim of 0.8 ± 1.9 kJ mol −1 originated from the inaccuracy of the interatomic potential functions employed for calculations of CED. An uncertainty of 1.9 kJ mol −1 for Δh mix corresponds to an expected error of the FH interaction parameter FH of about ± 3 in case of PEO-H 2 O mixtures at 300 K and a PEO weight fraction of 0.8. Such error exceeds values of the FH parameters for many common polymer-solvent combinations. [63] As a consequence, the parameterization of accurate thermodynamic models exclusively based on CED and the energy of mixing is virtually impossible, although the calculated SPs show a very good agreement with experiments (see Figure 2). Furthermore, the estimated error of Δh mix,sim rationalizes the contradicting results of solubility predictions using the FH theory for various polymer drug/active mixtures obtained recently, which only partially agree with experimental observation. [8,17,18] However, if specific interactions do not dominate the physicochemical properties of polymer mixtures, as indicated by relatively low values of el , the calculated contributions vdw and el to SP can be used for qualitative solubility predictions potentially facilitating rapid screening of the thermodynamic compatibility between polymers and biologically active compounds (see Figure 3).
Contrary to the FH theory, the parameterization of the PHSC EOS used not only the temperature dependence of CED but also of pressures. This enabled a more accurate modeling of the intermolecular interactions by employing an LJ potential using both, pure and mixed states of the compounds under consideration. The subsequent separation of the intermolecular potential into repulsive and attractive contributions facilitated the targeted parameterization of the reference part of the PHSC EOS Z ref = 1 + Z HS + Z CB representing an ensemble of noninteracting HS chains and the perturbation part Z pert including all attractive interactions. In this way, considerably improved thermodynamic modeling of aqueous PEO solutions was achieved compared to the FH theory providing fairly good agreement of Δh mix,EOS with experiments as well as negative Δs mix,EOS . The decrease of entropy due to mixing of two compounds showing strong hydrogen bonding is connected with the formation of complexes in solution as observed, e.g., for alcohol-water mixtures. [29] Similarly, previous theoretical studies [31] on the thermodynamic properties of PEO-H 2 O mixtures found negative Δs mix , which coincides largely with the results obtained in this work. In addition, negative Δs mix was, e.g., obtained for mixtures of meloxicam, a nonsteroidal antiinflammatory drug, and aqueous PEO solutions. [32] This emphasizes the need for accurate thermodynamic models including the effects of pronounced hydrogen bonding on the physicochemical properties of polymer-active mixtures. Even though the PHSC EOS parameterized in this work provided a reasonable thermodynamic description of PEO-H 2 O solutions, its accuracy is still insufficient to reproduce both, CED and pressures as a function of temperature, with equal quality. However, by using more sophisticated expressions for Z pert and applying more accurate models for consideration of specific interactions as discussed above, the present computational methodology appears a promising way for predictions of the capability of polymers to solubilize biologically active compounds. Moreover, the PHSC EOS has been extended for the thermodynamic modeling of statistical, block, and alternating copolymers as well as their corresponding mixtures. [44] Therefore, the combination of such EOS with the atomistic simulation procedures presented here would provide a promising tool for the efficient, in silico design of sophisticated (co-)polymeric nanocarriers.

Conclusions
In summary, computational approaches for predictions of the thermodynamic compatibility between polymers and low molecular weight compounds based on atomistic models have been assessed for applications in nanomedicine. The IS sampling has been evaluated for predictions of physicochemical properties of www.advancedsciencenews.com www.advtheorysimul.com amorphous polymers, in particular CED and pressures. Sufficient accuracy of the IS sampling has been confirmed by comparison with comprehensive test calculations using the DS of polymer conformations. Main advantage of the IS sampling over the DS is its trivial parallelization.
The CED obtained from atomistic simulations along with their energetic contributions of electrostatic and van der Waals interactions have been used for calculations of the SPs based on the concept of Hansen. Together with the FH theory, the SPs have been employed for qualitative solubility predictions of various polymer-solvent combinations in good agreement with experimental data. However, the FH theory failed for polymer mixtures showing pronounced specific interactions such as hydrogen bonding.
In contrast, the PHSC EOS parameterized using atomistic simulations yielded good agreement with experimental data for PEO-water and PLA-water mixtures, including the enthalpy and entropy of mixing. However, this can be partially attributed to error compensations since PHSC EOS is still too inaccurate to reproduce simulated CED and pressures with the same quality. Consequently, accurate and reproducible predictions of phase diagrams could not be obtained since only small errors of the free energy lead to large deviations of the calculated phase diagrams. Nonetheless, the PHSC EOS provides reproducible values of the free energy of mixing that are expected to be sufficiently accurate for reliable predictions of the thermodynamic compatibility between polymers and small molecules even for systems showing strong specific interactions. Moreover, the PHSC EOS can be straightforwardly extended to copolymers for modeling of more sophisticated polymeric materials.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.