Group contributions to the solvation free energy from MST continuum calculations

Group contributions to the solvation free energy from MST continuum calculations

Ignacio Soteras^I; Antonio Morreale^II; José María López^I; Modesto Orozco^II; F. Javier Luque^I

^IDepartament de Fisicoquímica, Facultat de Farmàcia, Universitat de Barcelona, Avgda Diagonal 643, Barcelona 08028, Spain

^IIMolecular Modelling and Bioinformatic Unit, Institut de Recerca Biomèdica, Parc Científic de Barcelona, Josep Samitier 1-5, Barcelona 08028, Spain, and Departament de Bioquímica i Biologia Molecular, Facultat de Química, Universitat de Barcelona, Martíi Franquès 1, Barcelona 08028, Spain

ABSTRACT

Group contributions to the free energy of solvation in water and octanol as well as to the octanol/water partition coefficient have been determined from Miertus-Scrocco-Tomasi continuum calculations. Particular attention is paid to the influence exerted by the procedure used to carry out the charge normalization in the MST model, as well as to the formalism of the partitioning scheme. A good agreement is found between the group contributions calculated by using different charge normalization and partitioning schemes in a series of structurally related drug-like molecule. Finally, the transferability of the group contributions determined for common chemical fragments along the series of molecules is discussed.

1 Introduction

The relevance of solvation in modulating the biological activity of drugs is well known for decades [1-4]. In fact, solvation influences the activity of a drug at two different-complementary levels. First, the drug must be soluble in aqueous environments, but it has to pass biological membranes which are strongly apolar. Accordingly, a proper balance between the global hydrophilic and lipophilic properties of the drug is required to ensure a correct distribution profile in the organism [5]. Second, binding to the target receptor implies that the drug (and the receptor) must be partially or totally desolvated. This energetically unfavourable process must be counterbalanced by the formation of suitable contacts between the drug and the receptor, which suggests that there must exist some complementarity between the electrostatic, steric and hydrophobic properties of the drug and the receptor binding site [6]. Therefore, not only the total free energy of solvation, but also the tridimensional distribution of polar and apolar regions along the molecular skeleton should be considered to understand biological activity of molecules [7].

The three dimensional representation of solvation can be obtained by means of the concept of group contribution [7]. These parameters are typically derived from empirical schemes. In some of these methods the contribution of a given fragment is calculated from the difference between the property values for related molecules in a homologous series, i. e. a compound bearing that fragment and the unsubstituted (reference) compound [8-10]. In the others, the molecules are subdivided into different fragments (atoms, functional groups) and their contributions are determined by using a suitable fitting procedure [for instance, see Refs. 11-13]. Those fragment-additive schemes assume that each fragment contribute a constant amount to the thermodynamic quantity, which can be estimated by adding the contributions of the different groups in the molecule. However, the near-independence on the molecular environment might also limit their use for the quantitative description of the solvation around molecules.

In the last years theoretical methods have been developed for the calculation of fragment contributions to the solvation free energy, particularly in the framework of quantum mechanical (QM) continuum solvation methods [14,15]. Thus, fractional methods based on the GB/SA methods have been developed by Cramer and Truhlar [16,17]. Our own group has developed also similar algorithms using the solvation method developed by Miertus-Scrocco-Tomasi (MST) model [18,19]. Though these methods demand a sizeable computational effort for large-sized molecules, they have the main advantage that the contributions are computed considering explicitly the molecular environment.

Following our previous studies, we examine here the influence of the MST-charge normalization procedure on the computed group contributions. Attention is also paid to the comparison between two different partitioning schemes defined in the context of NDDO-based QM MST calculations. Finally, the transferability of the group contributions within a series of structurally related drug-like molecules is also discussed.

2 Theoretical framework

2.1 The MST continuum model

The solvation process is divided in three stages [20]: i) generation of a cavity in the solvent, ii) insertion of the ''uncharged'' solute, and iii) generation of the solute charge distribution. Accordingly, the solvation free energy, DG_sol, can be computed as the addition of cavitation (DG_cav), van der Waals (DG_vW ) and electrostatic (DG_ele ) components (Eq. 1).

The cavitation free energy is determined following Pierotti's scaled particle theory [21] adapted to molecular shaped cavities by using the procedure proposed by Claverie [22]. Thus, the cavitation free energy of atom i, DG_{C - P,i}, is determined weighting the contribution of the isolated atom, DG_P,i, by the ratio between the solvent-exposed surface of such an atom, S_i, and the total surface of the molecule, S_T, as noted in Eq. 2, where N is the number of atoms.

The van der Waals term is computed using a linear relationship to the solvent-exposed surface of each atom, as noted in Eq 3, where DG_vW,i is the van der Waals free energy of atom i, and x_i is the atomic surface tension, which is determined by fitting to the experimental free energy of solvation [23-25].

Finally, the electrostatic contribution is determined assuming that the solvent is a continuum polarizable medium, which reacts against the solute charge distribution. Following the original formalism of the Pisa's Polarizable Continuum Model [14, 26], the reaction field generated by the solvent is introduced as a perturbation operator,

_R, into the Schrödinger equation (Eq. 4). Such a perturbation operator is expressed in terms of a set of imaginary charges located on the solute cavity (Eq. 5), which are obtained by solving the Laplace equation with suitable boundary conditions (Eq. 6).

where M is the total number of surface elements, j, in which the solute cavity is divided and q_j denotes the set of charges (located at r_j) that represents the solvent response.

where V_T is the total electrostatic potential, which includes both solute and solvent contributions, n is the unit vector normal to the surface element j, S_j is the area of the surface element j, and e is the solvent dielectric constant.

Finally, the electrostatic component of DG_sol is determined as noted in Eq. 7, where the index ''sol'' means that the perturbation operator is adapted to the fully relaxed charge distribution of the solute in solution, and the index ''o'' stands for the gas phase environment.

2.2 Charge normalization

A delicate issue in the MST (and PCM) model is the treatment of the charge normalization, which is necessary owing to two different effects: i) the solute charge density that lies outside the cavity, and ii) the numerical errors due to tessellation of the cavity surface. Accordingly, the sum of the apparent surface charges, Q_s = q_j , do not satisfy the relationship given by Eq. 8, where Q_M denotes the total charge of the solute, and the error in the charges is given by D_s = Q_theor - Q_s , where Q_theor = Q_M.

The original MST model [26] takes into account the charge normalization by using the expressions given in Eqs. 9a-b, where Q and Q denote the sum of positive and negative apparent charges, respectively.

Further efforts carried out by the Pisa group [27, 28] have given rise to a most elaborate treatment of the charge compensation, which allows us to treat separately the two sources of errors and presumably is better suited to account for local effects related to the asymmetry of the electron charge distribution. The new approach relies on two main features. First, the introduction of separate factors for the apparent charges induced by the nuclei (f^Z) and the electrons(f^e), which makes possible the direct correction of the numerical errors on the nuclear part due to the discretization of the cavity (Eq. 10). Second, the explicit calculation of the charge density outside the cavity (Q_out) which is determined from the flux of the solute electric field through the cavity surface (Eq. 11). By calculating with accuracy Q_out, and by taking into account its effect through an extra set of apparent charges, the solvent response induced by the electron charge density inside the cavity (Q_in) is then only affected by the numerical errors due to tessellation of the cavity (Eq. 12).

where E_M (s) is the electric field created at each tessera of the cavity.

2.3 Partitioning of the free energy of solvation

Partitioning of the non-electrostatic terms into atomic contributions is straigthforward, since they are related to the solvent-accessible surface of atoms see Eqs. 2 and 3). The partition of the electrostatic term is more difficult, but can be achieved by using a perturbation treatment of the mutual polarization between solute and solvent, which allows us to rewrite Eq. 7 as noted in Eq. 13 [29].

Since Eq. 13 computes the electrostatic component from the interactions between the charge distribution of the solute with the set of apparent charges that define the solvent reaction field (see Eq. 5), it allows us to decompose DG_ele into atomic contributions, as noted in Eq. 14. This procedure, which can be implemented in both semiempirical and ab initio versions of the MST model, computes the fractional contributions to DG_ele of a given atom ifrom the interaction energy between the whole charge distribution of the molecule with the apparent charges located at the surface elements pertaining to the portion of the cavity generated from that atom. This scheme will be denoted hereafter surface-based partitioning method.

where c stands for the basis of atomic orbitals, and P_mn is the mn element of the one-electron density matrix.

In the framework of NDDO-based methods, an alternative partitioning can also be envisaged, as noted in Eq. 15. Accordinhly, the fractional contribution of a given atom i is determined from the interaction energy between the elements mn associated to the atom i with the whole set of apparent charges spread over the cavity. This approach, which has already been examined by treating the solute charge distribution at the classical level [30], is denoted atom-based partitioning scheme.

Whatever the method chosen to partition the electrostatic component, the solvation free energy can be then expressed in terms of atomic contributions, DG_sol,i(Eq. 16).

2.4 Technical details of the simulations

Calculations were performed for a series of compounds closely related to tacrine (i.e., the first drug approved by the FDA for the treatment of Alzheimer's disease [31]; see Fig. 1). This series of molecules was chosen because it allows us to analyze the group contributions for three distinct fragments. First, the central 4-aminopyridine unit, which is common to all the compounds. Second, the alicyclic Z-(CH₂)₃ fragment, with Z being CH₂ or O (denoted by a and b, respectively, in the numbering of the compounds). Third, the aromatic X-Y-C(R')-C(R'') fragment, which supports all the chemical differences between compounds, mainly involving the presence of several heteroatoms. These chemical differences are exploited here to also examine the transferability of the contributions to the free energy of solvation and to the partitioning coefficient exerted by the common groups along the series of molecules.

MST calculations were performed in water and in octanol, which allows us to discuss the differential trends of the solvation in the two solvents. The standard parameters of the MST model [23-25] were used in all the calculations. Geometries of the compounds were optimized at the HF/6-31G(d) and AM1 levels in the gas phase and kept frozen in MST calculations. For the surface-based partitioning scheme, calculations were done at both HF/6-31G(d) and AM1 levels. With respect to the atom-based approach, calculations were carried out only at the AM1 level. MST calculations were carried out using local versions of MOPAC [32] and MonsterGauss [33] programs.

3 Results and discussion

The electrostatic component, DG_ele, of the free energy of solvation in water and octanol for the series of tacrine-like compounds determined from MST-HF/6-31G(d) and MST-AM1 calculations are shown in Table 1. As expected from the presence of polar groups, there is a significant contribution of the electrostatic term to the solvation of neutral compounds in water (DG_ele values varying from -9 to -20 kcal/mol). Compared to the results obtained in water, the electrostatic free energy in octanol is reduced by around 60%, reflecting the lower polarity of octanol compared to water.

The HF/6-31G(d) results determined by using the two charge normalization procedures are very similar, as noted in the small magnitude of the mean-signed error (mse < 0.4) and root-mean square deviation (rmsd < 0.3), as well as in the scaling coefficient of the regression equations, which is close to unity (see Table 1). In fact, the differences arising from the two charge normalization schemes are clearly smaller than those due to the use of ab initio and semiempirical versions of the MST continuum model, especially for the solvation in water. This finding is also clearly reflected in the electrostatic component of the octanol/water partition coefficient (log; see Eq. 17), as noted in the results shown in Table 2, which evidence the close similarity in the log values obtained from the two charge normalization procedures. In summary, no big changes can be expected in the total free energy values due to the use of different charge normalization procedures.

The group contributions to DG_ele and log of the three structural fragments are given in Table 3. Due to the chemical modifications in the aromatic X-Y-C(R')-C(R'') fragment, the contribution of this fragment is very different along the series of compounds (for instance, it varies from around -1 to -7 kcal/mol for the hydration free energy). In contrast, a roughly constant contribution is found for the alicyclic Z-(CH₂)₃ unit, with average values of -0.6 (-0.3) and -4.2 (-1.8) kcal/mol for the solvation in water (octanol) of compounds 1a-8a (Z=CH₂) and 1b-8b (Z=O), respectively. In turn, the contribution to log of the Z-(CH₂)₃ fragment is rather constant (see Table 3), it being on average 0.3 for Z=CH₂ and -1.6 for Z=O (in logP units). Finally, the average contribution of the 4-aminopyridine unit, which is common to all the molecules, to the solvation in water amounts to -8.4 kcal/mol, though it varies in a range of near 4 kcal/mol for the different compounds. A similar behaviour is found for the solvation in octanol (average contribution of -3.5 kcal/mol), where it varies in a range of near 3 kcal/mol. As a result, the 4-aminopyridine fragment has an average contribution of 3.6 logP units, though the values obtained for the series of compounds vary by around 1 logP unit. These findings reveal the large influence played by the chemical environment, specially for the hydration free energy, which makes it necessary to be cautious about the transferability of the group contributions to the solvation free energy.

Table 3

Based on the preceding results, it can be stated that the two charge normalization formalisms behave similarly for the series of tacrine-like compounds examined here, which gives support to the conclusions derived previously for selected small organic compounds [28]. Moreover, its influence on the group contributions is sensibly smaller than the effect exerted by the chemical environment. Nevertheless, the appearance of significant differences between the two procedures for molecules having a larger anisotropy in the charge distribution might not be ruled out.

Let us now turn our attention to the dependence of the group contributions on the scheme adopted to partition the solvation free energy. As noted above, two different partitioning schemes of DG_ele can be defined in the context of MST-AM1 calculations. In the surface-based approach (Eq. 14), the fractional contribution of a given atom is determined from the interaction of the whole molecule with the solvent's reaction field generated in the solvent-exposed surface of such an atom. On the contrary, in the atom-based approach (Eq. 15) it is determined from the interaction of the atom with the solvent's reaction field in the entire surface of the molecule.

In principle, for compounds having particular structural features the two partitioning schemes might provide very different results. However, in practice, for the series of tacrine-like compounds, our results indicate that the two partitioning schemes provide very similar group contributions to both DG_ele and log (see Table 4). Thus, the largest differences amount to around 0.7 and 0.5 kcal/mol for the fragmental contributions to the solvation in water and in octanol, and to around 0.3 logP units for the partition coefficient. This agreement strongly supports the qualitative value of the fractional contributions derived here, since they do not appear to be very dependent on the arbitrary selection of a given partitioning scheme [30].

To investigate the transferability of the group contributions to both DG_ele and log , each of these quantities (X) was calculated (see Eq. 18) by adding the contribution of the X-Y-C(R')-C(R'') fragment to the average contribution (X)determined for the 4-aminopyridine unit, which is common to all the molecules, and for the Z-(CH₂)₃ fragment, which is shared by compounds 1a-8a (Z=CH₂) and 1b-8b (Z=O). The thermodynamic quantities estimated from Eq. 18 were subsequently compared with the values determined for the whole molecule from QM MST-HF/6-31G(d) calculations.

The results of the statistical analysis are shown in Table 5. Overall, the mse between the fragment-additive and QM values is close to zero in all cases. Moreover, the fragment-additive and QM results correlate very well (0.91 < r < 0.95), and the scaling coefficients are close to unity. However, the rmsd amounts to 0.9 and 0.6 kcal/mol for the solvation in water and octanol, thus reflecting the nonnegligible effect due to the chemical environment on the fragmental contribution of the 4-aminopyidine unit (see above), which is the common structural unit along the series of compounds. This effect, nevertheless is less apparent for the electrostatic component of the octanol/water partition coefficient, since the rmsd between the fragment-additive and QM values is reduced to only 0.3 (in logP units).

4 Conclusion

Group contributions determined for the three basic units of the tacrine-like compounds examined here are reasonably similar irrespective of the charge normalization procedure and of the partitioning scheme used. In fact, the results point out that the dependence of the group contribution for the 4-aminopyridine unit (i.e. the common structural unit along the series of compounds) on these two factors is clearly smaller in magnitude than the effect played by the chemical environment. Though the procedure used to renormalize the apparent surface charge does not have major effect in the results for the tacrine-like molecules, caution is necessary since more relevant differences might appear for molecules with larger anisotropy in the charge distribution. Finally, the similarity found between the results determined from the two partitioning schemes gives us confidence in the use of this theoretical approach to explore the transferability of fragmental contributions to the solvation free energy and to the partitioning between different solvents. This information can be valuable to design with confidence parameters for QSAR studies of bioactive molecules.

Acknowledgements

We thank Prof. Tomasi for a version of his original PCM code, which was modified to implement our ab initio and semiempirical versions of the MST method. This work has been supported by the Ministerio de Ciencia y Tecnología (BIO2003-06848 and SAF2002-04282), and the Centre de Supercomputació de Catalunya.

Received on 22 September, 2003

• [1] M. Orozco and F. J. Luque, Chem. Rev. 100, 4187 (2000)
• [2] C. J. Cramer and D. G. Truhlar, Chem. Rev. 99, 2161 (1999).
• [3] W. L. Jorgensen, Acc. Chem. Res. 22, 184 (1989).
• [4] P. A. Kollman, Chem. Rev. 93, 2395 (1993).
• [5] C. A. Lipinski, F. Lombardo, B. W. Doming, and P. Feeney, Adv. Drug Deliv. Rev. 23, 3 (1997).
• [6] Y. C. Martin, Quantitative Drug Design. A Critical Introduction, Marcel Dekker, New York (1978).
• [7] C. Hansch and A. Leo, Exploring QSAR: Fundaments and Applications in Chemistry and Biology, American Chemical Society, Washington DC (1995).
• [8] A. Leo, C. Hansch and D. Elkins, Chem.Rev., 71, 525 (1971).
• [9] C. Hansch and A. Leo, Substituents Constants for Correlation Analysis in Chemistry and Biology, Wiley, New York (1979).
• [10] R. F. Rekker, The hydrophobic fragmental constant, Elsevier, New York (1977).
• [11] S. Cabani, P. Gianni, V. Mollica, and L. Lepori, J. Solut. Chem., 10, 563 (1981).
• [12] A. K. Ghose, V. N. Viswanadhan, and J. Wendoloski, J. Phys. Chem. A, 102, 3762 (1998).
• [13] R. Wang, Y. Fu, and L. Lai, J. Chem. Inf. Comput. Sci., 37, 615 (1997).
• [14] J. Tomasi and M. Persico, Chem. Rev. 84, 2027 (1994).
• [15] J. L. Rivail and D. Rinaldi, In Computational Chemistry, Review of Current Trends, J. Leszcynski (Ed.), World Scientific Publishing, Singapore (1996), pp 139-169.
• [16] C. J. Cramer and D. G Truhlar, Chem. Phys. Lett. 74, 198 (1992).
• [17] G. D. Hawkins, C. J. Cramer, and D. G. Truhlar, J. Phys. Chem. B 102, 3257 (1998).
• [18] F. J. Luque, X. Barril, and M. Orozco, J. Comput.-Aided Mol. Des. 13, 139 (1999).
• [19] J. Muńoz, X. Barril, B. Hernández, M. Orozco, and F. J. Luque, J. Comput. Chem. 23, 554 (2002).
• [20] F. J. Luque, C. Curutchet, J. Muńoz, A. Bidon-Chanal, I. Soteras, A. Morreale, J. L. Gelpí, and M. Orozco, Phys. Chem. Chem. Phys. 5, 3827 (2003)
• [21] R. A. Pierotti, Chem. Rev. 76, 717 (1976).
• [22] P. Claverie. In Intermolecular Interactions: From Diatomics to Biomolecules, B. Pullman (Ed.), Wiley, Chichester (1978).
• [23] F. J. Luque, Y. Zhang, C. Aleman, M. Bachs, J. Gao, and M. Orozco, J. Phys. Chem. 100, 4269 (1996).
• [24] F. J. Luque, C. Alemán, M. Bachs, and M. Orozco, J. Comput. Chem. 17, 806 (1996).
• [25] C. Curutchet, M. Orozco, and F. J. Luque, J. Comput. Chem. 22, 1180 (2001).
• [26] S. Miertus and J. Tomasi, Chem. Phys. 65, 239 (1982).
• [27] B. Mennucci and J. Tomasi, J. Chem. Phys. 106, 5151 (1997).
• [28] M. Cossi, B. Mennucci, J. Pitarch, and J. Tomasi, J. Comput. Chem. 8, 833 (1998).
• [29] F. J. Luque, J. M. Bofill, and M. Orozco, J. Chem. Phys. 103, 10183 (1995).
• [30] A. Morreale, J. L. Gelpí, F. J. Luque, and M. Orozco, J. Comput. Chem., 24, 1610 (2003).
• [31] X. Barril, S. G. Kalko, M. Orozco, and F. J. Luque, Mini-Rev. Med. Chem. 2, 27 (2002).
• [32] MOPAC 6.0. Version locally modified by F.J. Luque and M. Orozco. University of Barcelona. 2001.
• [33] M. Peterson and R. Poirier, MonsterGauss, Department of Biochemistry, Univ. Toronto, Canada, Version modified by R. Cammi and J. Tomasi (1987), and by F. J. Luque and M.Orozco (2002).

Publication Dates

• Publication in this collection
  27 Apr 2004
• Date of issue
  Mar 2004

History

• Received
  22 Sept 2003