Solubility calculation of pharmaceutical compounds – A priori parameter estimation using quantum-chemistry

Fluid Phase Equilibria 299 (2010) 161–170

Contents lists available at ScienceDirect

Fluid Phase Equilibria journal homepage: www.elsevier.com/locate/fluid

Solubility calculation of pharmaceutical compounds – A priori parameter estimation using quantum-chemistry Jan Cassens a , Feelly Ruether a,∗ , Kai Leonhard b , Gabriele Sadowski a a b

Laboratory of Thermodynamics, Department of Biochemical and Chemical Engineering, Technische Universitaet Dortmund, Emil-Figge Str. 70, D-44227 Dortmund, Germany Laboratory of Technical Thermodynamics, RWTH-Aachen, Schinkelstr. 8, D-52062 Aachen, Germany

a r t i c l e

i n f o

Article history: Received 18 August 2010 Received in revised form 17 September 2010 Accepted 20 September 2010 Available online 25 September 2010 Keywords: Solubility Pharmaceuticals PCP-SAFT Equation of state Quantum chemistry

a b s t r a c t Pure-component electrostatic properties for pharmaceutical compounds and intermediates (xanthene, ibuprofen, aspirin, p-hydroxyphenylacetic acid, p-toluic acid and o-anisic acid) were obtained by quantum-chemical methods. Afterwards, these properties were used for the a priori determination of the pure-component parameters for the Perturbed-Chain Polar Statistical-Associating Fluid Theory (PCP-SAFT). These parameters were applied to perform solubility calculations for binary solute–solvent mixtures. In these calculations the only parameter fitted was the binary parameter. The results show a good agreement of the modeled solubility and experimental data for the considered solutes in nonpolar and polar solvents. Finally, the application of different combination rules to also predict the binary interaction parameter in the mixture was investigated. © 2010 Elsevier B.V. All rights reserved.

1. Introduction The knowledge of the solubility of pharmaceuticals in pure solvents and solvent mixtures is crucial for designing separation or purification steps, especially crystallization processes. The choice of a good solvent usually depends on the amount of the solubility data available and requires a high experimental effort which is expensive and time consuming. Reducing this effort to a minimum is only possible by applying thermodynamic models which are applicable to calculate or even predict the solubility of a drug in solvents and even solvent mixtures. On the one hand, there exist empirical correlations, such as quantitative structure–properties relationships (QSPR), as proposed among others by Ghasemi and Saaidpour [1]. These empirical relationships use structure-dependent descriptors which are related to the aqueous solubility of the pharmaceutical. Ran and Yalkowsky [2] proposed the general solubility equation (GSE) for organic nonelectrolytes in water utilizing the water–octanol partition coefficient and the melting point as parameters for the solubility calculations. On the other hand, if the activity coefficient of the pharmaceutical solute is known, the solubility can be obtained from thermodynamic relationships. For this purpose, models like NRTL-SAC [3], COSMO-RS [4], and PC-SAFT [5] have already been

applied. The big advantage of the latter thermodynamic models is that – compared to the empirical methods – a solubility estimation can be provided based on much less experimental information. The applicability of the PC-SAFT equation of state for solubility calculations has already been shown by Ruether and Sadowski [6]. This model requires pure-component parameters for each component (solvents and solutes). These parameters are usually determined by fitting to pure-component vapor pressure and liquid density data. They can subsequently be used to calculate mixture properties by applying an additional binary parameter. Since for pharmaceutical solutes in most cases vapor pressures and liquid density data are not available, these parameters usually have to be fitted to binary (solubility) data as done in [6]. The aim of this work is to reduce the number of fitted model parameters and to determine at least some of them from quantum-chemical methods. 2. Thermodynamic model The solubility of a pharmaceutical can be obtained by a simplified thermodynamic phase–equilibrium relationship according to



xiL ∗ Corresponding author. Tel.: +49 231 755 3199; fax: +49 231 755 2572. E-mail address: [email protected] (F. Ruether). 0378-3812/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2010.09.025

SL −H0i 1 = exp i RT



1−

T SL T0i



(1)

Here xiL is the mole fraction in the liquid phase (L), T is the system temperature, R is the ideal gas constant and  i is the activity coef-

162

J. Cassens et al. / Fluid Phase Equilibria 299 (2010) 161–170

Nomenclature List of symbols a Helmholtz free energy per number of particles (J) (6) Cij dispersion coefficient ClP fij SL H0i k kij m Neff,i P qi Q ri R RC S T SL T0i V xi

regression constant conversion factor enthalpy of fusion (J/mol) Boltzmann constant 1.38065 × 10−23 (J/K) binary interaction parameter segment number effective number of electrons fitted parameter in regression charge of atom i (C) ˚ quadrupole tensor (DA) coordinate vector of atom i ideal gas constant 8.314 (J/molK) ˚ mean radius of curvature (A) exposed surface (A˚ 2 ) system temperature (K) melting temperature (K) volume (A˚ 3 ) mole fraction of component i

Greek letters ˛ shape factor of a spherical molecule ˛ ¯i static segment polarizability i activity coefficient of the component i ε/k dispersion energy parameter (K) εassoc /k association energy parameter (K) i scalar dipole moment of component i (D) assoc association volume parameter descriptor in regression ∗0 04∗ descriptor in regression ˚ ␪ traceless quadrupole moment (DA) ˚ i scalar qudrupole moment of component i (DA)  i,xx ,  i,yy ,  i,zz diagonal entries of the traceless quadrupole ˚ tensor (DA) ˚  segment diameter (A) Superscripts L liquid phase res residual hc hard-chain disp dispersion dipole dipole quadpol quadrupole assoc association Subscripts fit fitted i, j components i and j 0i pure substance i 0 first estimation

SL and T SL are the enthalpy ficient of the solute component i. H0i 0i of fusion and the melting temperature, respectively. These are pure-component properties and can be determined, e.g., by DSC measurements. For simplification and due to the lack of experimental data the influence of the heat capacity on the solubility in Eq. (1) is neglected. In this work, the activity coefficient is obtained using the PCPSAFT equation of state (EOS) [5,7,8]. Within this model, the residual

Helmholtz energy ares is calculated based on a perturbation expansion with the hard-chain fluid as reference according to ahc ares adisp adipol aquadpol = + + + kT kT kT kT kT

(2)

where ahc is the hard-chain contribution of m tangentially connected hard spheres, each of them having a segment diameter . adisp accounts for the dispersive interactions between the chain segments [5]. Helmholtz-energy contributions due to multipolar interactions are accounted for by adipol for dipole–dipole interactions and by aquadpol for quadrupole–quadrupole interactions as proposed by Gross and Vrabec [7,8]. To describe a pure substance, the PCP-SAFT EOS according to Eq. (2) requires five pure-component parameters: besides the parameters m and  these are the dispersive interaction energy parameter ε as well as the scalar dipole moment  and the scalar quadrupole moment . The last-mentioned two substance properties, accounting for the polar interactions, can be taken either from literature or can be derived from quantum-chemical calculations as shown earlier by Leonhard et al. [9,10]. Using an additional contribution to the Helmholtz energy, associative effects can also be accounted for by an additional term aassoc [11], which leads to: ares ahc adisp adipol aquadpol aassoc = + + + + kT kT kT kT kT kT

(3)

For the parameterization of the associative interaction, the association-energy parameter εassoc as well as the association volume assoc have to be defined for each associating component [11]. These parameters have also to be fitted for solvents and solutes to experimental vapor-pressure and solubility data, respectively. Within this work, the dipole moments and the quadrupole moments as well as the static polarizabilities for the solvent components have been taken from literature. If these properties were not available in literature, they have been calculated with the Gaussian03 software package [12]. Electrostatic properties of the solute components were all calculated with the Gaussian03 software package [12]. Since the tensors given by Gaussian03 are in the standard-orientation of the molecule, they were transformed into the tensor’s principle axis frame. Moreover, the center of the coordinate system had to be transformed from the center of nuclear charges, as used in Gaussian03, to the center of mass which is more common for EOS applications. The transformation of a quadrupole moment tensor with a non-vanishing dipole moment is described by Gray and Gubbins [13]. Furthermore, in EOS the traceless quadrupole moment  is more commonly used. It is defined as [14] =

3 2



Q−



1 Tr(Q )1 3

=

1 qi (3r i r i − r 2i 1) 2

(4)

i

In Eq. (4) Tr(Q) is the trace of the quadrupole tensor Q and terms the sum over all elements on the main diagonal. In this context traceless means that the sum over all diagonal elements of Q is zero. The symbol qi is the charge of atom i with coordinates ri . In the PCP-SAFT EOS the multipole moments are used as scalar values. For this reason the length of the dipole vector was used for dipolar interactions. The scalar quadrupole moment  i was calculated from the quadrupole moment tensor of molecule i in its principal axis 2 , 2 2 frame having non-vanishing diagonal elements i,xx and i,zz i,yy [9]: 2 2 2 2 + i,yy + i,zz ) = i2 ( 3 i,xx

(5)

J. Cassens et al. / Fluid Phase Equilibria 299 (2010) 161–170

163

2.1. Estimation of electrostatic properties For the calculation of the electrostatic properties, possible conformers of each component were identified and a geometrically optimized configuration was generated. For these calculations, a TZVP basis set and the B3LYP method were used. Afterwards, the electrostatic properties were calculated with an aug-cc-pVDZ basis set and the second-order Møller–Plessert method. For a detailed description of the basis sets and the methods used for the calculation the reader is referred to the Gaussian03 manual [12]. As mentioned before, the quadrupole – given in the standard orientation – has been transferred to the principal axis frame with the center of mass as origin. As different states of energy have been identified for the considered conformers, a Boltzmann average has been calculated over the different energetic states to average the electrostatic properties. The geometrical properties were directly derived from the fused hard-sphere model based on the geometrically optimized structure. The dipole dispersion coefficients can be calculated by using frequency-dependent polarizabilities [10]. In this work the Slater–Kirkwood approximation [15] was used for the estimation of the C(6) dispersion coefficients, where the dispersion coefficient between two atoms is defined as: (6) Cij

˛ ¯ i˛ ¯j 3 = 1/2 1/2 2 (˛ ¯ i /Neff,i ) + (˛ ¯ i /Neff,j )

(6)

The effective number of electrons Neff,i of component i depends on the static polarizability ˛ ¯ i and the dipole dispersion coefficient (6) Cii . (6)

Neff,i =

16 Cii 9 ˛ ¯i

(7)

If one assumes that the static polarizability as well as the effective number of electrons are additive properties, the dispersion coefficient of one component can be calculated by summing up all atom–atom dispersion coefficients. The effective number of electrons for several atoms and atomic dispersion coefficients have been developed by Wu and Yang [16] for the atoms H, C, O, and N. Dispersion coefficients for molecules up to C8 have been calculated by Wu and Yang who give the mean deviation of these values with 1%. The dispersion coefficients between two different atoms were estimated by inserting Eq. (7) into Eq. (6) which leads to the combination rule (6) (6) 2

(6) Cij

=

2 ] 2[(Cii Cjj ) Neff,i Neff,j (6)

2 ) (Cii Neff,i

1/3

(6)

1/3

2 ) + (Cij Neff,i

1/3

(8)

2.2. A priori estimation of pure-component PCP-SAFT parameters To estimate the PCP-SAFT pure-component size and shape parameters of solutes from quantum chemistry, a fused hardsphere model with element-specific radii was used to describe the atomic configuration of each substance as proposed by Leonhard et al. [17]. Knowing the chemical composition of a solute, an equilibrated structure of atoms was obtained by applying quantum-chemical calculations using the Gaussian03 software package. Based on this equilibrated structure, a fused hard-sphere model was constructed with each sphere representing one atom placed at the atoms coordinates. The geometrical properties of this model are well defined since the position of the spheres is given by the atom coordinates and the radii of the spheres correspond to the element-specific radii that were estimated in [17] as follows. To calculate the volume, the exposed surface and the shape factor of the fused hard-sphere model of the solute, a methodology proposed by Alejandre et al. [18] was used and afterwards these properties

Fig. 1. Schematic drawing of the molecule methane as a fused hard-sphere model. The radius of curvature is determined by averaging the radii from the center to a new supposed surface (left). The molecular volume and the exposed surface (right) are calculated from the exposed part of the molecule.

were related to the PCP-SAFT parameters m and . For a certain training set of molecules, Van Nhu et al. [17] fitted the PCP-SAFT pure-component parameters to experimental vapor–liquid equilibrium and density data. These fitted parameters will be referred to hereafter as mfit and  fit . The volume of the fused hard-sphere model corresponds to the volume Vfit calculated from the PCP-SAFT pure-component parameters mfit and  fit . It is defined as [17] Vfit =



fit 4 mfit 3 2

3

(9)

The shape factor ˛ of the fused hard-sphere model describes the deviation of the molecule from a spherical structure. It is calculated from the mean radius of curvature RC of the molecule, the exposed surface S and the molecular volume V which are in turn obtained from the fused hard-sphere model [18]: ˛=

RC S 3V

(10)

The shape factor of the fused hard-sphere model can be related to the shape factor calculated with the fitted segment number of the PCP-SAFT EOS by [19] ˛fit =

mfit + 1 2

(11)

Fig. 1 shows a schematic drawing of the fused hard-sphere model for methane and illustrates the derivation of the geometrical properties. Based on the fused hard-sphere model, a new surface is generated like it is shown on the left side of Fig. 1. For this step the molecule is wrapped in a foil-like surface. The mean radius of curvature is the average over all radii Ri from the center of the molecule to the assumed surface. The grey surface on the right side of Fig. 1 indicates the exposed part of the fused hard-sphere model which is considered for the volume and the surface calculation. The volume Vfit and the shape factor ˛fit were determined for the components of the training set from PCP-SAFT parameters by Van Nhu et al. [17] according to Eqs. (9) and (11). Afterwards, the element-specific radii of the spheres of the fused hard-sphere model have been fitted iteratively to these properties by minimizing the deviation between the properties calculated by PCP-SAFT and those calculated by the fused hard-sphere model. The solvent molecules of the training set were composed of the atoms C, H, ˚ and O and the element-specific radii of these atoms (rH = 1.078 A, ˚ and rO = 1.453 A) ˚ were used in the calculations concernrC = 1.43 A, ing the solute molecules in this work. Thus, applying this approach a determination of the geometrical properties was possible for solute molecules consisting of C, H, and O atoms. From the fused hard-sphere model for each solute the volume V0 , the surface S0 , and the mean radius of curvature RC,0 were calculated with the element specific radii determined by Van Nhu et al. [17], respectively. The index “0” denotes the first-estimation values directly deduced from the fused hard-sphere model properties, which means that these properties were calculated based on the molecule geometry and the given element-specific radii. Using

164

J. Cassens et al. / Fluid Phase Equilibria 299 (2010) 161–170

these geometrical properties, first-estimation values for the segment diameter  0 , segment number m0 and dispersion energy ε0 as well as the reduced dipole moment ∗0 and the reduced quadrupole moment 04∗ , were determined according to Eqs. (12)–(16).

expressed by [20]

m0 = 2˛0 − 1

where fij is a conversion factor,  ij is the segment diameter



0 = 2

ε0 = ∗0 = 04∗ =

3

(12)

Cij

(20)

ij6 (6)

3 V0 4 m0

(14)

calculated according to Eq. (19) and Cij is the mixed dispersion coefficient per segment. Leonhard et al. [10] proposed an expression for mixed dispersion coefficients based on London’s approximation and the static segment polarizabilities ˛ ¯ i of the components i:

(15)

Cij ≈ 2

(13)

C (6) m20 06

(6) (6)

 m0 ε0 03 4

(16)

m0 ε0 05

A further discussion of the development and the relevance of each descriptor can be found in Ref. [17]. The final PCP-SAFT parameters P (standing for the PCP-SAFT parameters , m, or ε) are obtained from the first-estimation values by subsequently fitting to the descriptors compensating for the non-ideal character of the polar and the dispersive interactions. For this purpose a multi-linear regression has been used [17]: P P P P = ε0 CεP0 + m0 Cm + 0 CP0 + 04∗ C P4∗ + ∗0 C ∗ + C0 0 

0

To describe mixtures, combination rules are applied to the EOS parameters. For the dispersion energy parameter and the segment diameter, Lorentz–Berthelot combination rules were used in the original PC-SAFT model [5]:



εij = (1 − kij )

ii + jj 2

εii εjj

(18) (19)

In Eqs. (18) and (19)  ij is the segment diameter of the fluid mixture of components i and j, εij is the dispersion energy parameter between component i and j and kij describes the binary interaction parameter. Leonhard at al. [9] proposed an approach to predict the dispersion energy parameter in the mixture. Based on the assumption that the model for the dispersion contribution of the Helmholtz energy truly and only covers the dispersive interaction, the dispersion energy parameter between two spherical molecules εij can be

(6)

(6)

˛ ¯ 2j Cii + ˛ ¯ 2i Cjj

˛ ¯ i˛ ¯j

(21)

The approach is similar to an earlier one proposed by Kohler for unsegmented molecules [21]. Inserting Eq. (20) in Eq. (21) and assuming the conversion factors to be molecule-independent constants, the London formula results in: εij ≈

2

εii ii6 εjj jj6 ˛ ¯ i˛ ¯j

(22)

ij6 εii ii6 ˛ ¯ 2j + εjj jj6 ˛ ¯ 2i

Another expression for the combined dispersion energy parameter can be deduced by substituting Eq. (20) in Eq. (18) which leads to εij ≈

2.3. Combination rules

Cii Cjj

(6)

(17)

0

In Eq. (17), Ci P mean the constants for the calculation of parameter P and descriptor i, respectively. To determine these constants, the first-estimation parameters and descriptors have been calculated for a training set of molecules from the fused hard-sphere model by Van Nhu et al. [17]. Afterwards they fitted the constants of Eq. (17) by inserting the first-estimation parameters in Eq. (17) and minimizing the deviation of the regressed parameter P to the one fitted to pure-component data. In this work we used the regression constants determined by Van Nhu et al. [17]. The regression in Eq. (17) is necessary because neither the PCPSAFT model nor the descriptors are exact. E.g., the presently used model neglects the dipole–quadrupole and the polar induction interactions, which is compensated by fitting the dispersion energy parameter. Such interactions are not accounted for in the first calculation of the dispersion energy parameter ε0 from the fused hard-sphere model. Therefore, the final pure-component parameters are optimized according to Eq. (17).

ij =

(6)

εij = fij

ii3 jj3 ij6

(6)

f

ij

fii fjj

Cij





(6) (6)

Cii Cjj

εii εjj

(23)

Under the assumption that the mixed conversion factor is the geometric mean of the pure-component conversion factors, which is a weaker assumption for the factors than being identical, Eq. (23) becomes [9] εij ≈

ii3 jj3 ij6

(6)



Cij

(6) (6)

Cii Cjj



εii εjj

(24)

Eqs. (22) and (24) offer two possibilities to replace the binary interaction parameter by explicit mathematical expressions. In Eq. (22) the static polarizabilities are required to define the dispersion energy parameter of the mixture. In Eq. (24) the dispersion coefficients are required to determine the dispersion energy parameter of the mixture. The static polarizabilities can be taken from quantumchemical calculations directly whereas the dispersion coefficients have to be calculated using dynamic polarizabilities [10] or estimated, as described in Section 2.1. 3. Results 3.1. A priori estimation of pure-component parameters The estimation of the PCP-SAFT pure-component parameters was performed for the solutes xanthene, ibuprofen, p-toluic acid, phydroxyphenylacetic acid, o-anisic acid, and aspirin, respectively. After the electrostatic properties have been calculated for the pharmaceutical molecules, the pure-component PCP-SAFT parameters for those molecules were a priori calculated as described above. The results are listed in Table 1. The physical properties like melting temperature and enthalpy of fusion of the solutes, as they are given in Table 1, have been taken from literature. The pure-component parameters for the solvents used in this work were fitted to experimental vapor–liquid equilibrium data [25]. If available, values for electrostatic properties and dispersion coefficients were also taken from literature. The PCP-SAFT parameters for the solvents considered in this work are given in Table 2.

J. Cassens et al. / Fluid Phase Equilibria 299 (2010) 161–170

165

Table 1 A priori calculated pure-component solute parameters. Component

m (−)

Aspirin Ibuprofen o-Anisic acid p-Hydroxyphenylacetic acid p-Toluic acid Xanthene a b c d

3.247 4.166 2.998 3.151 2.983 3.186

˚  (A)

ε/k (K)

4.085 3.998 4.029 3.981 3.964 4.196

 (D)

294.821 249.984 290.765 268.32 281.142 345.969

3.557 1.7 2.237 2.653 2.946 1.102

˚  (DA)

˛( ¯ A˚ 3 )

25.695 11.724 12.959 8.248 14.623 13.139

C(6) (a.u.) d

18.104 24.696 15.851 15.956 15.698 22.805

5439.89 7238.71d 4221.93d 4006.51d 3797.33d 7670.71d

SL

H0i (J/mol) a

25600 21900b 19710 a 28000 c 22719a 19200a

SL T0i (K)

408.15a 349.15b 374.15a 755.85c 452.75a 373.65a

Data from literature [22]. Data from literature [23]. Data from literature [24]. Values calculated with Eq. (8).

Table 2 Pure-component parameters for the solvents. Component Acetone Benzene Cyclohexane Acetic acid Ethanol Methanol Toluene Water a b c d

˚  (A)

m (−) a

ε/k (K) a

2.529 2.008a 2.486a 1.757 4.546 3.361 2.460a 1.910

3.390 3.957a 3.855a 3.712 2.5542 2.528 3.920a 2.417

˚  (DA)

 (D) a

236.962 303.647a 281.383a 315.999 193.8257 219.095 300.190a 306.559

a

2.912 0a 0a 2.012 1.441b 1.698b 0.374a 1.861b

a

4.525 7.992a 0.729a 9.7568 5.749b 4.274b 7.481a 2.895b

˛( ¯ A˚ 3 )

C(6) (a.u.) a

786.660 1772.650a 2269.340a 555.003d 521.036d 218.216d 2557.720a 45.290c

6.420b 10.215b 10.485b 5.069 5.110b 3.270b 12.167 1.458b

Values are taken from literature [17]. Values are taken from literature [10]. Values are taken from literature [26]. Values calculated with Eq. (8).

If no reference is given for the electrostatic properties in the table, these values have been calculated with Gaussian03. It has to be noted here, that for ethanol and methanol high values for the segment number were determined compared to the parameters already published in literature [11] where these substances were modeled as associating fluids. However, in this work the a priori estimation of pure-component parameters was performed neglecting the associative contribution in the PCP-SAFT EOS since this effect was not considered in the development of the regression equations according to Eq. (17). Hence, in the first part of the solubility calculations in Section 3.4, the pure-component parameters as listed in Tables 1 and 2 as well as Eq. (2) were utilized in the modeling. In the second part, the associative contribution was taken into account (Eq. (3)) and additional association parameters for the solute molecules as presented in Table 3 were used for the solubility calculations. 3.2. Determination of association parameters As some of the molecules in this study exhibit associative interactions, this effect should be taken into account in the modeling [6,27,28]. Moreover, solvents like ethanol and water are usually modeled as associating components [11]. However, for the estimation of the solute PCP-SAFT parameters the association was Table 3 Pure-component association parameters for associating solutes. Component

εassoc /k (K)

Aspirin 2197.634 Ibuprofen 811.377 o-Anisic 1243.797 acid p627.517 Hydroxyphenylacetic acid p-Toluic 1455.851 acid

assoc (−)

Proton donors/proton acceptors

0.0001 0.1659 0.002

1/1 1/1 1/1

0.0269

2/2

0.0013

1/1

initially neglected since there is presently no reliable predictive method to calculate the association parameters. To include the association, it is assumed that the a priori parameter calculation covers the dispersive interaction as truly dispersive. Thus, the purecomponent parameters for the dispersive and polar interaction of the solutes were taken from the quantum-chemical calculations (Table 1). The two additional association parameters for the solute molecules εassoc and assoc have been fitted subsequently to the experimental solubility data of binary solute–solvent systems. Table 3 shows the fitted pure-component association parameters for the associating solutes. For each considered solute molecule one proton donor and one proton acceptor site has been assumed for the calculations except for p-hydroxyphenylacetic acid where two proton donor and two proton acceptor sites were assumed. The pure-component parameters for the associating solvents ethanol, methanol and water were taken from Gross and Sadowski [11]. The association parameters for solvents are in the range of 103 –104 K for the association energy parameter and 10−4 –10−1 for the association volume parameter, respectively. The values of the association parameters are in the same ranges as reported in literature for other associating molecules [11]. The difference in the quality of the calculated solubility using the approach according to Eq. (2) compared to the approach according to Eq. (3) will be presented and discussed in Section 3.4. 3.3. A priori combination rules Leonhard et al. [9] proposed an approach for a priori combination rules and applied it to binary systems of small (solvent) molecules. They used the combination rules according to Eqs. (22) and (24) to calculate the dispersion energy parameter of the mixture. This approach was also applied to the pharmaceutical systems considered in this work. For the calculations we used the purecomponent properties as given in Tables 1 and 2. The results of Eq. (22) – referred to as London combination rule – have been attained using the segment-based static polarizabilities. Applying the dispersion energy parameter of the mixture εij from

166

J. Cassens et al. / Fluid Phase Equilibria 299 (2010) 161–170

Table 4 Comparison of a priori calculated and fitted binary interaction parameters for the solutes xanthene, p-hydroxyphenylacetic acid (p-hpa) and ibuprofen in selected solvents. Component Xanthene p-Hpa Ibuprofen

Solvents

C(6)

London

kij -fitted

Benzene Cyclohexane Toluene Toluene

0.01778 −0.00956 −0.00422 −0.17231

0.00754 0.01982 0.00071 0.02265

0.01218 0.02537 0.06934 0.01087

Eq. (22), the corresponding binary parameters kij , as defined by the Lorentz–Berthelot mixing rule in Eq. (18), can be estimated. These binary interaction parameters for three solutes in selected solvents are compared in Table 4 to the ones calculated by Eq. (24) – referred to as C(6) combination rule – and additionally to those from the fitting results to solubility data – referred to as kij -fitted. For the fitting of the binary interaction parameters, one data point (the last or the first) of the respective experimental dataset has been used. The solvents considered for the calculations are non associating solvents which is in accordance with the assumption that the dispersive interaction is supposed to be truly dispersive and not influenced by other interactions. For the solutes p-hydroxyphenylacetic acid and ibuprofen the influence of self-association has been neglected in these cases and the parameters according to Table 1 have been used. The results in general show a smaller deviation of the fitted binary parameters to those calculated by London’s combination rule compared to the ones calculated with the C(6) combination rule. This shows that the values calculated by London’s combination rule could be used as a rough estimation of the binary interaction parameter for the systems given in Table 4. The changing sign of the binary interaction parameters obtained from the C(6) combination rule can be ascribed to inaccuracies of using the Slater–Kirkwood approximation for the calculation of the dispersion coefficients. This effect possibly occurs due to the assumption that static polarizability as well as the effective number of electrons are additive properties and have been estimated with a group contribution method [16]. The influence of the various interaction parameters on the quality of solubility calculations will be shown in Section 3.5. The binary interaction parameters have also been calculated for the systems containing the solvents water, ethanol, and methanol according to Eqs. (22) and (24). For these calculations the purecomponent parameters according to Tables 1 and 2 have been

utilized. The binary interaction parameters were calculated with a high deviation (0.13–0.3) compared to the fitted ones. These significant deviations indicate that not only dispersive but other interactions (e.g., association) which are not accounted for explicitly in Eq. (2) are relevant for those systems and thus have to be compensated by the binary interaction parameter. Thus, it has to be concluded that the a priori combination rules do not work for the systems including water, ethanol or methanol. For this reason the results for these systems are not presented in Table 4. 3.4. Solubility calculations using the a priori pure-component parameters The a priori estimated pure-component parameters for the pharmaceutical components were used to model their solubility in various solvents. The binary interaction parameters were fitted here to each binary system. Fitting of this parameter was performed to only one solubility data point (the first or the last data point in the experimental data series) as well as to two solubility data points (the first and the last point in the experimental data series). The average deviation of modeled and measured solubility could not always be improved significantly by fitting to more than one data point, as can bee seen in Table 5. On the other hand, this means that only one data point was sufficient to determine the binary parameter and to obtain a satisfying description of the experimental data. Table 5 presents the average relative as well as the average absolute deviation of calculated and experimental solubility for the considered solutes. For the solubility of ibuprofen in toluene and in acetone the average relative deviation is only 1.95% and 0.5%, respectively. The solubility of xanthene and p-hydroxyphenylacetic acid in several solvents can also be modeled in good agreement with experimental data. However, as expected from the previous chapters, the calculated solubility for the solutes p-toluic acid, o-anisic acid, and aspirin in water deviates significantly from experimental data. This is again assigned to the fact that the association between the water molecules as well as between water and solute molecules is not explicitly taken into account in the PCP-SAFT EOS using the a priori determined parameters. Figs. 2–4 exemplarily show solubility diagrams for the solutes xanthene and ibuprofen in different solvents for the considered temperature range. In the diagrams the experimental data are depicted as symbols whereas the lines represent the modeled solubility. It can be seen that the temperature dependence of the solubility can be well described by the PCP-SAFT model using the

Table 5 Average deviation of simulation results from the experimental values for several systems. The binary interaction parameters have been fitted to one data-point and to two data-points, respectively. System

Average relative deviation (kij fitted to one data-point) (%)

Average relative deviation (kij fitted to two data-points) (%)

Xanthene/benzenea Xanthene/cyclohexanea Ibuprofen/tolueneb Ibuprofen/methanolb Ibuprofen/ethanolb Ibuprofen/acetoneb p-Toluic acid/water a p-Hydroxyphenylacetic acid/toluenec p-Hydroxyphenylacetic acid/waterc o-Anisic acid/watera Aspirin/watera

5.86 5.47 1.91 22.99 9.51 0.5 >30 4.38

5.7 5.4 2.1 22.4 9.2 0.5 >30 4.3

10.49

8.3

a b c

>30 >30

Experimental data from literature [22]. Experimental data from literature [23]. Experimental data from literature [24].

>30 >30

Average absolute deviation (mole fraction)

Temperature range (K)

Number of datapoints

Experimental solubility (mole fraction)

2 × 10−2 3 × 10−2 4 × 10−3 4 × 10−2 2 × 10−2 1 × 10−3 6 × 10−5 3 × 10−4

299.6–349.6 331.7–359.7 283.15–308.15 283.15–308.15 283.15–308.15 283.15–308.15 288.35–356.45 283.15–288.15

6 6 5 5 5 5 16 2

0.16–0.61 0.21–0.72 0.1–0.3 0.1–0.3 0.12–0.33 0.14–0.32 3.6 × 10−5 –4.3 × 10−4 1.03 × 10−3 –1.25 × 10−3

6 × 10−5

283.15–298.15

4

3.49 × 10−3 –7.14 × 10−3

2 × 10−4 4 × 10−4

278.15–346.15 278.15–345.15

20 15

1.86 × 10−4 –6.21 × 10−3 1.99 × 10−4 –1.21 × 10−3

J. Cassens et al. / Fluid Phase Equilibria 299 (2010) 161–170

Fig. 2. Solubility (mole fraction over temperature) of xanthene in cyclohexane (kij = 0.02537) and benzene (kij = 0.01218). The pure-component parameters for the solute were calculated a priori. Experimental data [22] are given by symbols and simulation results are represented by lines.

Fig. 3. Solubility (mole fraction over temperature) of ibuprofen in toluene (kij = 0.01087) and acetone (kij = −0.01376). The pure-component parameters for the solute were calculated a priori. Experimental data [23] are given by symbols and simulation results are represented by lines.

a priori pure-component parameters as well as only one constant binary parameter for each solute–solvent system. Fig. 5 shows the comparison of calculated and experimental solubility of aspirin in water. Also for this calculation the PCP-SAFT parameters as listed in Tables 1 and 2 were used. The modeling underestimates the solubility in the lower temperature range. However, it has to be noted here that for this system the absolute values of solubility are very small as they are in the range of 10−4 . The absolute error is also in the range of 10−4 .

167

Fig. 4. Solubility (mole fraction over temperature) of ibuprofen in ethanol (kij = 0.02734). The pure-component parameters for the solute were calculated a priori. Experimental data [23] are given by symbols and simulation results are represented by lines.

Fig. 5. Solubility (mole fraction over temperature) of aspirin in water (kij = −0.03362). The pure-component parameters for the solute were calculated a priori. Experimental data [22] are given by symbols and simulation results are represented by lines.

In general it can be concluded that the method used for parameter estimation can be utilized to characterize solute substances consisting of C, H, and O-atoms within the PCP-SAFT framework. The so-obtained parameters can be applied to calculate the solubility of the considered solutes in different non-associating solvents. The deviation of the calculated solubilities from the experimental data is in an acceptable range for most systems. However, for systems containing ethanol, methanol, and especially for aqueous systems significantly higher deviations from the experimental data are observed.

Table 6 Average deviation of simulation results from the experimental values for considered systems. The approach according to Eq. (2) is compared to the approach including associative interactions according to Eq. (3). Binary systems

Average relative deviation using Eq. (2) (%)

Average relative deviation using Eq. (3) (%)

Ibuprofen/methanol Ibuprofen/ethanol p-Toluic acid/water p-Hydroxyphenylacetic acid/water o-Anisic acid/water Aspirin/water

23.0 9.5 >30 10.5 >30 >30

4.91 4.4 28.4 4.7 16.7 >30

168

J. Cassens et al. / Fluid Phase Equilibria 299 (2010) 161–170

Fig. 6. Solubility (mole fraction over temperature) of ibuprofen in methanol. Comparison of results accounting for the association term in the model (kij = 0.00177, dotted line) and simulation results modeled according to Eq. (2) (kij = −0.00519, solid line) to experimental data [23] (circles).

3.5. Solubility calculations including the associative interactions To investigate the influence of the associative contribution as described in Section 3.2, two additional association parameters had to be fitted subsequently for the solvent components. Hence, at least two additional experimental data points for the fitting of the association parameters εassoc and assoc were required. The fitting was done to the available experimental solubility data of the binary systems. For the solvent components the parameters were taken from Gross and Sadowski [11]. Due to the application of new solvent parameter sets, the binary interaction parameters had to be refitted for each binary system. Table 6 compares the results for the modeling using these new parameters accounting for association with those used earlier neglecting the association (see also Table 5). For ibuprofen in the solvents methanol and ethanol as well as for o-anisic acid in water the accuracy of the simulated solubility is significantly improved. For other aqueous systems, only a slightly lower average relative deviation of the calculated solubility from the experimental data could be observed. As an example, Fig. 6 compares the calculated solubility of ibuprofen in methanol with and without accounting for association. The calculated solubility using the association term is in a good agreement to the experimental data. It can be seen that the result can be remarkably improved by taking into account the association contribution explicitly in the modeling. The results for the solubility of o-anisic acid in water are depicted in Fig. 7. The comparison of the modeling results to the experimental data shows that in this case considering the associative effects explicitly can only slightly improve the quality of the modeling in the lower temperature region. Concerning the solubility of aspirin in water no significant change in modeling accuracy can be observed by including the associative interactions in the model. Fig. 8 shows the comparison of the modeling results with the approach with association and the one without associative contribution to the experimental data. In both cases the slope of the simulated solubility curve is too steep at higher temperatures whereas the solubility at lower temperatures is underestimated. In conclusion the comparison between the approaches (accounting for association and neglecting the associative contribution, respectively) shows that for the solute ibuprofen in the solvents methanol and ethanol considering the associating interaction improves the modeling results significantly. For aqueous

Fig. 7. Solubility (mole fraction over temperature) of o-anisic acid in water. Comparison of results accounting for the association term in the model (kij = −0.05485, dotted line) and simulation results modeled according to Eq. (2) (kij = −0.04310, solid line) to experimental data [23] (circles).

systems only a slight or even no improvement can be observed. Moreover, it has to be noted that the inclusion of the association term requires two additional model parameters for each component which have to be fitted to experimental data. 3.6. Pure predictions of solubility Figs. 9 and 10 show a comparison of the predictive results obtained using the a priori calculated pure-component parameters from Tables 1 and 2 as well as the a priori calculated binary parameters from Table 4. Hence, no experimental data point was required to attain the results presented in Figs. 9 and 10. The binary systems xanthene in benzene as well as ibuprofen in toluene, respectively, were chosen as representative examples. As it can be seen for the system xanthene in benzene, the a priori calculated pure-component parameters in combination with the binary parameters obtained with London’s combination rule or C(6) combination rule can be successfully used to predict the xanthene solubility. However, for the solubility of ibuprofen in toluene, the C(6) combination rule clearly overestimates the experimental data. London’s combination rule underestimates the solubility but its deviation is much smaller than the one of the C(6) combination rule. The results for the latter system show that the a priori combination rules are often but not in general capable to predict the binary system with sufficient accuracy. A reasonable explanation

Fig. 8. Solubility (mole fraction over temperature) of aspirin in water. Comparison of results accounting for the association term in the model (kij = −0.17246, dotted line) and simulation results modeled according to Eq. (2) (kij = −0.03362, solid line) to experimental data [23] (circles).

J. Cassens et al. / Fluid Phase Equilibria 299 (2010) 161–170

Fig. 9. Solubility (mole fraction over temperature) of xanthene in benzene. Comparison of results with fitted kij (gray line), with kij calculated using London rule (black line) and with kij calculated using C(6) rule (black dashed line) to experimental data [22] (symbols).

Fig. 10. Solubility (mole fraction over temperature) of ibuprofen in toluene. Comparison of results with fitted kij (gray line), with kij calculated using London rule (black line) and with kij calculated using C(6) rule (black dashed line) to experimental data [23] (symbols).

for the deviation between the calculated and the fitted binary mixing parameters could be the fact that the predictive combination rules presented here were developed to describe pure dispersive interactions between the molecules only and the assumption was made, that all other contributions (e.g., polar contributions) are correctly covered by separate terms. 4. Conclusion In the first part of the work the solubility of several pharmaceutical solute components was modeled in different solvents using the PCP-SAFT EOS where the dispersive and the polar interactions were taken into account without consideration of the associative contribution. For each solute the dipole moment, the quadrupole moment as well as the static polarizability were calculated with the Gaussian03 software package. The dispersion coefficients were estimated with the Slater–Kirkwood approximation. Using quantum-chemical descriptors which were estimated from the electrostatic and geometric properties, a priori pure-component parameters for the PCP-SAFT EOS were derived. These parameters were capable to model the solubility of the solutes xanthene, p-hydroxyphenylacetic acid and ibuprofen in non-associating sol-

169

vents in fair agreement with experimental data. For systems containing methanol and ethanol as well as for aqueous systems the relative deviation of calculated and experimental data was higher than 10%. However, the absolute deviation for the aqueous systems lies in the range of 10−4 which causes a high relative error at a small absolute deviation. In general, this approach provides an elegant alternative to fitting all pure-component parameters to binary solubility data. It provides a well-defined parameter set which is usually difficult to obtain from fitting procedures. In the second part, the association was explicitly accounted for in the solubility calculation of systems containing associating solute and solvent components, especially water. By doing so the quality of the solubility calculations for the aqueous system aspirin–water did not change remarkably whereas the modeling of the systems ibuprofen–methanol and ibuprofen–ethanol could be improved significantly. However, taking association into account requires at least two more experimental solubility data points for fitting the pure-component parameters of the associating component. Finally, it was shown that the London a priori combination rule can be used as a rough estimation for binary interaction parameters for modeling the solubility in the systems xanthene in benzene and cyclohexane as well as for ibuprofen and p-hydroxyphenylacetic acid in toluene. In contrast the application of the binary interaction parameter calculated by the C(6) combination rule results in higher deviations from experimental values. However, both combination rules do not always provide reliable values of the binary interaction parameters for solubility calculations. Consequently, the binary interaction parameter between solute and solvent should better be determined by fitting to experimental data. In this case, already fitting to one data point only leads to satisfactory modeling results which cannot be remarkably improved by more extensive fitting. Eventually, the proposed method provides the opportunity for the a priori estimation of pure-component solute parameters for the PCP-SAFT model. Only one binary interaction parameter has to be fitted to one experimental data point, hence the experimental effort can be reduced significantly. Compared to the common approaches of parameter fitting, it delivers well-defined parameter sets on a reduced data-basis. Acknowledgment Part of this work was performed within the Cluster of Excellence “Tailor-Made Fuels from Biomass”, which is funded by the Excellence Initiative by the German federal and state governments to promote science and research at German universities. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]

J. Ghasemi, S. Saaidpour, Chem. Pharm. Bull. 55 (2007) 669–674. Y.Q. Ran, S.H. Yalkowsky, J. Chem. Inf. Comput. Sci. 41 (2001) 354–357. C.C. Chen, P.A. Crafts, Ind. Eng. Chem. Res. 45 (2006) 4816–4824. A. Klamt, F. Eckert, M. Hornig, M.E. Beck, T. Burger, J. Comput. Chem. 23 (2002) 275–281. J. Gross, G. Sadowski, Ind. Eng. Chem. Res. 40 (2001) 1244–1260. F. Ruether, G. Sadowski, J. Pharm. Sci. 98 (2009) 4205–4215. J. Gross, J. Vrabec, AIChE J. 52 (2006) 1194–1204. J. Gross, AIChE J. 51 (2005) 2556–2568. K. Leonhard, N. Van Nhu, K. Lucas, Fluid Phase Equilib. 258 (2007) 41–50. M. Singh, K. Leonhard, K. Lucas, Fluid Phase Equilib. 258 (2007) 16–28. J. Gross, G. Sadowski, Ind. Eng. Chem. Res. 41 (2002) 5510–5515. M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, J.A. Montgomery Jr., T. Vreven, K.N. Kudin, J.C. Burant, J.M. Millam, S.S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G.A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J.E. Knox, H.P. Hratchian, J.B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, P.Y. Ayala, K. Morokuma, G.A. Voth, P. Salvador, J.J. Dannenberg, V.G. Zakrzewski, S. Dapprich, A.D. Daniels, M.C. Strain, O. Farkas, D.K. Malick, A.D. Rabuck, K. Raghavachari, B. Foresman, J.V. Ortiz, Q. Cui, A.G. Baboul, S. Clifford, J. Cioslowski, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-

170

[13] [14] [15] [16] [17] [18] [19] [20]

J. Cassens et al. / Fluid Phase Equilibria 299 (2010) 161–170 Laham, C.Y. Peng, A. Nanayakkara, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, M.W. Wong, C. Gonzalez, J.A. Pople, Gaussian03, Revision C.02, Gaussian Inc., Wallingford, CT, 2004. G.G. Gray, K.E. Gubbins, Theory of Molecular Fluids, Clarendon Press, Oxford, 1984. K. Lucas, Molecular Models for Fluids, Cambridge University Press, Cambridge, 2007. T.A. Halgren, J. Am. Chem. Soc. 114 (1992) 7827–7843. Q. Wu, W.T. Yang, J. Chem. Phys. 116 (2002) 515–524. N. Van Nhu, M. Singh, K. Leonhard, J. Phys. Chem. B. 112 (2008) 5693–5701. J. Alejandre, S.E. Martinezcasas, G.A. Chapela, Mol. Phys. 65 (1988) 1185–1193. J.M. Walsh, K.E. Gubbins, J. Phys. Chem. 94 (1990) 5115–5120. K. Lucas, Angewandte Statistische Thermodynamik, Springer-Verlag, Berlin/Heidelberg, 1986.

[21] F. Kohler, Chem. Monthly 88 (1957) 857–877. [22] J. Marrero, J. Abildskov, Solubility and Related Properties of Large Complex Chemicals, DECHEMA, Frankfurt/Main, 2003. [23] J. Abildskov, Organic Solutes Ranging from C2 to C41, DECHEMA, Frankfurt/Main, 2005. [24] S. Gracin, A.C. Rasmuson, J. Chem. Eng. Data 47 (2002) 1379–1383. [25] T.E. Daubert, R.P. Danner, Data Compilation Tables of Properties of Pure Compounds, American Institute of Chemical Engineers, New York, 1985. [26] S.M. Cybulski, T.P. Haley, J. Chem. Phys. 121 (2004) 7711–7716. [27] I. Tsivintzelis, I.G. Economou, G.M. Kontogeorgis, J. Phys. Chem. B 113 (2009) 6446–6458. [28] G.K. Folas, S.O. Derawi, M.L. Michelsen, E.H. Stenby, G.M. Kontogeorgis, Fluid Phase Equilib. 228 (2005) 121–126.