Automatic parametrization of non-polar implicit solvent models for the blind prediction of solvation free energies

Solvation, in which separated solvent and solute molecules are combined to form a solution, is an elementary physical process in nature that provides a foundation for more complicated chemical and biological processes, such as ion channel permeation, protein-ligand binding, electron transfer, signal transduction, DNA specification, transcription, post-transcription modification, gene expression, and protein synthesis. Therefore, the understanding of solvation is a prerequisite for the quantitative study of other more complex chemical and biological processes and is of paramount importance to chemistry, physics, and biology.^14,67,78 The most basic and reliable experimental observation of solvation is the solvation free energy, which measures the free energy change in the solvation process. As a result, one of the most important tasks in the solvation modeling and computation is the prediction of solvation free energies, which has recently drawn much attention in computational biophysics and chemistry.^{14,15,46,67,78,88}

A large variety of solvation models has been developed in the past few decades for solvation analysis and solvation free energy prediction. In general, these models can be classified into either knowledge based models or physics based models. Knowledge based models usually utilize a large amount of available data of solvation free energies to train statistical models for the solvation free energy prediction. One example of this type of models makes use of solvent accessible surface areas, solvent excluded surface areas, or calibrated atomic surface areas for solvation free energy predictions.⁸⁵ Knowledge based models can provide relatively accurate predictions provided that enough experimental data are available. On the other hand, physics based solvation models are based on fundamental laws of physics. These models can be further classified into three major classes. One of them is explicit solvent models in which both solvent and solute molecules are described at the atomic or electronic level. This class of models is accurate in general, but is usually computationally expensive.^45,49,72 Another class of models is integral equation based solvation theories,^36,65,66 in which the solute molecule is still modeled at the atomic level, while the solvent is modeled by statistical mechanics, such as liquid density functional theory, Ornstein-Zernike equation with hypernetted-chain, Percus-Yevick equation, or other specified closure. Integral equation based solvation models reduce the number of degrees of freedom dramatically compared to explicit solvent models. A major feature of integral equation approaches is that these methods are able to provide a good approximation of solvent microstructures near the solvent-solute interface.²⁸ The other class of the physics based solvation models is implicit solvent models,^{1,7,12,22,27,61,64,71,89} in which the solute molecule can be modeled at either the atomic or the quantum mechanical level, while the solvent is described as a dielectric continuum.

Compared to other methods, implicit solvent models have advantages of involving a small number of degrees of freedom, highly accurate modeling of strong and long-range electrostatic interactions that often dominate solvation phenomena, and convenient treatments of solvent-solute electronic polarization effects.¹² In implicit solvent models, the solvation free energy is typically divided into polar and non-polar components. The polar solvation free energy is due to electrostatic interactions and can be modeled by a number of approaches, including Born dielectric sphere model, Kirkwood model,⁴³ generalized Born (GB) model,^{2,16,19,30,31,41} polarizable continuum model (PCM),^18,39,81 and Poisson-Boltzmann (PB) theory.^{29,37,62,70,73} Among these approaches, Born dielectric sphere model which treats the solute molecule as the dielectric sphere with a centrally located atomic charge is the simplest one. The Kirkwood model gives analytical expressions for electrostatic potentials in a spherical solute-solvent system with multiple charges inside the solute molecule.⁴³ The GB model^{2,16,19,30,31,41} is able to deal with arbitrarily shaped molecules and account for the impact of each point charge by its effective distance (i.e., the Born radius) from the solvent-solute interface.^13,80 The GB model is relatively faster, while depending on other methods, such as the PB theory, for its parametrization. The PCM describes solute electrons by using quantum mechanics so that solvation induced polarizations can be accounted for in an iterative procedure. A few different versions of PCM models have been reported, including dielectric PCM, integral equation formalism PCM, and the variational PCM models.^18,39,81 PCM is relatively computationally expensive due to its quantum mechanical charge calculation while its accuracy is bounded by the PB type of treatment of electrostatic potential as well as the quality of solvent-solute interfaces. The PB theory can be derived from more fundamental Maxwell’s equations.^{27,37,62,70,73} It is more accurate than the GB model and computationally more efficient than the PCM. The PB can be easily coupled with the quantum mechanics for a more accurate description of the solute charge density.^9,86,87 It has become one of the most popular solvation models due to its relatively low computational cost and high modeling accuracy for biomolecules.

The non-polar solvation component has been modeled by a number of terms. A popular approach, based on the scaled-particle theory (SPT) for non-polar solutes in aqueous solutions,^63,76 is to use a solvent-accessible surface area (SASA) term.^53,79 It was shown that a solvent-accessible volume (SAV) term is relevant in large length scale regimes.^38,51 Recent studies indicate that SASA based solvation models may not describe van der Waals (vdW) interactions near solvent-solute interface.^11,22,23,82 A combination of surface area, surface enclosed volume, and vdW potential has been shown to provide accurate non-polar solvation predictions.^10,83

The polar and non-polar components are typically decoupled in traditional implicit solvent models. Recently, new efforts have been made to couple polar and non-polar components in the literature.^7,17,89 Differential geometry based solvation model makes use of the differential geometry theory of surfaces to dynamically couple polar and nonploar models by surface evolution, which is driven by free energy optimization.^7–9,89 Coupled with an optimization procedure, this model was shown to deliver some of the best non-polar solvent free energy analysis.¹⁰ By applying constrained optimization to non-polar parameter selection, this model provides state-of-the-art solvation free energy fitting and cross validation results for a large number of solute molecules.⁸³ 

It is important to validate solvation models by experimental data. SAMPLx blind solvation prediction project, aimed at testing protein ligand binding models, also provides benchmark tests of the performance of the solvation models. Many different solvation approaches have been tested by the SAMPLx challenging molecules.^{21,25,26,34,35,55,56,58–60,68,69} Implicit solvent models are among the most competitive approaches in SAMPLx tests.

The objective of the present work is to develop an automatic protocol to parametrize implicit solvent models for the blind solvation free energy prediction. In our approach, we utilize PB theories for computing electrostatic solvation free energies. Both a point charge based PB model and a Kohn-Sham density functional theory (KSDFT) based polarizable PB model are developed in the present work. For the KSDFT based PB model, an iterative process is developed to take care of the solvent polarization and solute response.⁹ In our non-polar models, we employ either a collection of surface area, molecular volume and van der Waals interactions, or a combination of atomic surface areas and molecular volume. To determine parameters in our non-polar models, we assume that similar molecules admit the same set of non-polar parameters. In this work, we develop two parametrization algorithms. One of them is based on a functional group scoring, and the other is based on a nearest neighbor search. Both algorithms utilize machine learning methods to optimize parameters. A database contains the solvation free energies in aqueous solvent for 668 solute molecules which is collected from the literature^4,57 to validate the proposed parametrization protocol.

The rest of this paper is organized as follows. In Section II, we present our solvation models. Both polar and non-polar models are discussed. Sections II B and II C present two non-polar parametrization algorithms, namely, functional group scoring and nearest neighbor approaches. An algorithm for unified blind solvation free energy prediction is developed in Section II D. The blind prediction of a large number of solvation free energies is carried out in Section III. The leave-one-out test of 668 molecules is employed to validate the proposed protocol. Finally, the accuracy and robustness of the proposed protocol is demonstrated by five SAMPL test sets. This paper ends with a conclusion.

As an implicit approach, our solvation model has polar and non-polar components,

where ΔG is the total solvation free energy, ΔG^p is the polar or electrostatic solvation free energy, and G^np is the non-polar solvation free energy. The related polar and non-polar solvation models are described in Secs. II A 1 and II A 2.

For polar solvation modeling, we utilize either the standard Poisson model or a KSDFT based Poisson model (i.e., a polarizable Poisson model) for polar solvation free energy calculation. These two approaches are described below.

The Poisson model is given by Refs. 27 and 73,

where ϵ(r) is the discontinuous dielectric constant and is chosen as ϵ_m = 1 in the solute domain Ω_m and at ϵ_s = 80 in the solvent domain Ω_s. Here ϕ is the electrostatic potential, N_m is the total number of partial charges in the solute, and q_i is the ith charge at the position of r_i. Additionally, ϕ_s and ϕ_m are limiting values of the potential from solvent and solute domains, respectively. Here n is the outer norm vector at point r of the interface Γ. Finally, Ω and ∂Ω are computational domain and domain boundary, respectively.

We use the solvent excluded surface (SES) in the present work and it is generated by using our in-house Eulerian solvent excluded surface (ESES) software.⁵⁰ ESES automatically embeds the SES into the Cartesian mesh, which is convenient for the finite difference type of Poisson solvers. The matched interface and boundary based PB (MIBPB) software package^6,27,84 is used to solve Eq. (2) with appropriate charge force fields for {q_i} and radius sets in SESs.

The polar solvation free energy is then calculated as

where ϕ_home(r_i) is the solution of the Poisson equation with uniform dielectric constant ϵ_m over the whole domain Ω.

Instead of using a given partial charge force field, the KSDFT based polarizable Poisson model utilizes a DFT approach to compute the charge distribution in the solvent environment. As such, this approach partially takes into account for solvent polarization and solute response in an iterative procedure. The combination of the Poisson equation with SES and KSDFT has been reported in the literature.^86,87 Recently, we have proposed a differential geometry based polarizable PB model that further allows solvent-solute interface to be optimized in the Poisson-Kohn-Sham (PKS) formalism.⁹ In the present work, we adopt SES in our PKS software for the sake of an accurate estimation of atomic surface areas via our ESES software package.⁵⁰ Consequently, our MIBPB software is also employed for solving the Poisson equation with sharp interfaces.^6,27,84

In this approach, the Poisson model is⁹ 

where ρ_total(r) = qn(r) − qn_n(r) is the charge distribution in the solute domain, with q being the unit charge of an electron, n(r) being the electron density, and $n n ( r ) = ∑ I Z I δ ( r − R I )$ being the nucleus density, in which Z_I and R_I are the atomic number and position vector of nucleus I, respectively.

The Poisson equation (4) is coupled to the following generalized KS equation:⁹ 

where m is electron mass and ħ² = h/2π with h being the Planck constant. Here ψ_i and E_i are the ith KS orbital and ith eigenvalue, respectively, and $U eff ≐q ϕ RF + U eff 0$ is the effective potential of the generalized Kohn-Sham equation with ϕ_RF = ϕ − ϕ₀ being the reaction field potential of the solute in the solvent environment, in which ϕ and ϕ₀ are the electrostatic potentials generated by the solute in the solvent and solute environments, respectively. Finally, $U eff 0$ is the effective Kohn-Sham potential in the vacuum environment. Note that the solution of Eq. (5) gives rise to the electron density $n ( r ) = ∑ i | ψ i | 2$⁠, which, in turn, updates the charge distribution for Eq. (4). A self-consistent iterative process has been developed to solve Eqs. (4) and (5) to their convergence⁹ starting from solving the KSDFT in vacuum as an initial guess to the coupled system.

For KSDFT calculations, we employ the SIESTA solver,⁷⁵ which uses pseudopotential to eliminate the complicated effects of core electrons. In all calculations, the default double-ζ plus single polarization (DZP) basis is used. The MeshCutOff is set as 125 Rydberg and local density approximation (LDA) is used to approximate the exchange correlation potential. The solution method is set to be “diagon.”

In this formalism, the polar solvation free energy is computed by

For the non-polar solvation component, we use either vdW term based model or an atomic surface area based model.

In our earlier work, the non-polar solvation free energy was modeled by^7–9,89

where γ and A are, respectively, surface tension and area of the solute molecule. Additionally, p and V are, respectively, surface enclosed volume and hydrodynamic pressure difference. Finally, ρ₀ is the solvent bulk density, and U^vdW is the vdW interaction potential, i.e., the Lennard-Jones (LJ) potential, the detailed implementation of vdW interaction potential can be found in the work,^7,83

where N is the total number of atomic types in a given solute molecule and σ_s is the probe radius. Here, ε_α and σ_α are the vdW potential well depth and vdW radius of the αth type of atoms, respectively, and r_i is position of the ith atom in the α type of atoms.

In Eq. (7) the integration is over solvent domain Ω_s. This non-polar solvation free energy model was shown to offer excellent predictions of experimental data for various non-polar molecules.¹⁰ In our recent work, we have demonstrated superb results for a large number of polar and non-polar molecules when the non-polar model in Eq. (7) is combined with a PB model for the polar solvation component.⁸³ 

In this non-polar solvation free energy model, the parameters γ, p, σ_s, and ε_α need to be learned from a training set. However, one can set ρ₀ = 1 because its effect can be absorbed to ε_α.

The LJ potential offers a physical model for dealing with vdW interactions near the solvent-solute interface. However, a drawback of this non-polar term is that the LJ potential energy is nonlinear with respect to the probe radius (σ_s), which gives rise to a difficulty for a linear optimization process. Note that the non-polar energy is linear with respect to pressure and surface tension parameters.

Mathematically, for a small probe radius, the vdW term for a given atom is proportional to the atomic surface area. Therefore, atomic surface area approach used by Kollman and co-workers⁸⁵ should have a similar modeling effect as that of the LJ potential. Based on this observation, we propose the following non-polar solvation energy model:

where N is the total number of atomic types in a given solute molecule. Here, γ^α and Area^α are the surface tension and surface area of the αth type of atoms, respectively, and $Area i α$ is the surface area of the ith atom in the α type of atoms. In this non-polar solvation free energy model, the parameters γ_α and p need to be learned from a training set. Both the LJ and the atomic-surface-area based non-polar solvation models are validated in the present work.

In the above non-polar solvation models, the solvent excluded surface area, atomic solvent excluded surface areas, and molecular enclosed volume are calculated by our in-house software ESES.⁵⁰ 

In our non-polar solvation model, we assume that each functional group of molecules admits the same set of optimal parameters. Similar approach has been successfully used in the literature,⁷⁴ including ours.^10,83 In this work, we further incorporate machine learning type of methods for the non-polar solvation free energy prediction. To this end, the whole dataset, which contains 668 molecules, is partitioned into training sets and testing sets. In the leave-one-out test case, in each step we only leave one molecule out as the test set, and all the other molecules are regarded as the training set. In a specific leave-SAMPLx-out study, all the SAMPLx molecules are left as the testing set, while the remaining molecules are treated as the training set. Furthermore, the training set in this scenario is also divided into two parts: (i) the mono-functional group molecules, in which molecules of a specific functional group are used to train a set of parameters for all similar molecules and (ii) the poly-functional groups molecules in which scoring weights are used to weight the contributions of parameters of various involved functional groups.

Table I lists twenty functional groups used in the training set. We denote {T₁, T₂, …, T_n} as a given set of n molecules with a specific functional group from the above 20 groups. For a molecule indexed j, the solvation free energy calculated from the polar solvation free energy model, atomic surface areas and molecular volume are listed as

where j = 1, 2, …, n.

In our approach, the solvation free energy is modeled as the sum of the polar and non-polar parts, thus for the jth molecule, the corresponding modeled solvation free energy can be expressed as

where P≐{γ¹, γ², …, γ^N, p} is the set of parameters to be optimized.

For a given set of molecules with the same mono-func tional group, we denoted $ΔG ( P ) ≐ G 1 ( P ) , G 2 ( P ) , … , G n ( P )$⁠, and the associated experimental solvation free energy is denoted as $Δ G Exp ≐ Δ G 1 Exp , Δ G 2 Exp , … , Δ G n Exp$⁠. Then the optimal parameter set can be obtained by solving the following Tikhonov regularized least square problem (also known as ridge regression) which has a closed form solution:

where ‖∗‖₂ is the L₂ norm of the quantity ∗. Here λ is the regularization parameter chosen to be a large number, such as 100, in the present work to ensure the dominance of the first term and avoid over-fitting through controlling the magnitude of ‖P‖₂.

Most molecules in the SAMPL blind test sets involve poly-functional groups. In this case, we further employ the poly-functional group molecules in the training set for training the optimal relative scoring weights between different functional groups. According to the relative scoring weights, the scoring weights between all the functional groups can be obtained through a simple normalization procedure. We denote $T ̃ 1 , T ̃ 2 , … , T ̃ n ′$ as a given set of poly-functional group molecules that has the same functional groups in the training set. The associated optimized parameter sets are P₁, P₂, …, P_m, where m is the number of functional groups in this set, and each P_i, i = 1, 2, …, m is learned through solving the regularized optimization problem given by Eq. (12). For the jth molecule in this poly-functional group set, we model its solvation free energy as

where $‖ω ‖ 1 ≐ ∑ i = 1 m ω i =1$⁠, with ω_i being the scoring weight of ith functional group.

The relative scoring weights for the m functional groups associated to this set of poly-functional groups molecules can be learned via solving the following constraint optimization problem,

subject to

and

where $Δ G ̃ ( ω )$ and $Δ G ̃ Exp$ represent, respectively, the predicted and experimental solvation free energies for this group of poly-functional group molecules. Since both the optimization object given by Eq. (14) and the constraint conditions in Eqs. (15) and (16) are convex with respect to the scoring weights ω, the above constrained optimization can be easily solved via a convex optimization solver in the CVX software package.^32,33

In the rest of this section, we provide an example to illustrate the procedure of functional group based approach for solvation free energy prediction. As shown in Fig. 1, the target molecule 2-chlorosyringaldehyde (Pubchem ID: 53479) contains four different functional groups. We need to find scoring weights for these functional groups. Note that in the above molecular searching scheme, pairwisely, the relative weights for each two functional groups can be determined by solving the constrained optimization problem in Eqs. (14) and (15) for molecules in the two corresponding functional groups. According to the pairwise relative weights, the functional group scoring in the target molecule can be achieved.

The above functional group based model is tested to be able to provide very accurate solvation free energy prediction provided that for each test molecule, parameters for all of the involved mono-functional groups can be determined, and a set of molecules with the same ploy-functional groups can be found in the training set as well. However, for some complex molecules that contain functional groups beyond those listed in Table I, this method fails. For instance, considering a SAMPL1 test molecule, thifensulfuron, as shown in Fig. 2, it has a very complex structure and contains functional groups cannot be found in Table I. We therefore propose the following ranking algorithm.

To deal with general and complex molecules, we propose another approach to rank molecules and develop a nearest neighbor approach for non-polar solvation free energy prediction. Our molecular ranking is based on atomic features. Our basic assumption is that if two solute molecules have similar atomic features, their chemical properties are also similar. Therefore, these molecules should share the same parameter set for the non-polar solvation free energy model. This approach is further described below.

Our method is based on the ansatz that molecules chemical properties are mainly determined by atomic features, including structural features and atomic electrostatic features. Among them, atomic structural features include atomic type, atom hybridization state, and bonding information. Atomic electrostatic features include atomic charge, atomic dipole, atomic quadrupole, and atomic electrostatic solvation free energy. Atomic features are selected based the criteria that a given atomic feature is retained if it can effectively discriminate the previously mentioned mono-functional groups listed in Table I. Table II lists all the atomic features that are used for molecule ranking.

We also construct atomic features by using the statistics variables, i.e., average, maximum, minimum, and variation of quantities in Table II. Finally, we rule out redundant features, where redundancy is due to the fact that some atomic features are 100% correlated with each other. In this case, only one of these highly correlated features is retained.

Atomic features are calculated by following methods.

• Molecular structural information is parsed by the Open Babel software.⁵⁴ 

• Atomic charge, atomic dipole, and atomic quadrupole are obtained via the distributed multipole analysis (DMA) method,⁷⁷ in which the charge density is originally computed by the density function theory with B3LYP and 6-31G basis in the Gaussian quantum chemistry software.^3,20,47

• Atomic electrostatic solvation energy is calculated by our in-house MIBPB software.^6,27,84

With the above selected atomic features, intuitively, we measure the similarity of molecules by the Pearson correlation coefficient of atomic feature vectors. Specifically, we denote the features of two molecules as vector X≐{x₁, x₂, …, x_k} and vector Y≐{y₁, y₂, …, y_k}, then their similarity is measured by

where $x ̄ ≐ 1 k ∑ i = 1 k x i$ and $y ̄ ≐ 1 k ∑ i = 1 k y i$⁠, and k is the dimension of the feature space. The higher the correlation between two atomic feature vectors indicates the more similarity between molecules.

By using the above nearest neighbor based nearest molecule searching procedure, for a given molecule, we obtain the similarity ranking of the remaining molecules with this specified one, as well as correlations between molecules. In the nearest neighbor prediction of the non-polar solvation free energy, we learn the parameters in the non-polar solvation free energy expression, Eq. (9), by a given number of nearest molecules found in the training set. The number of the nearest molecules used is based on the following principles:

• All the molecules whose correlation with a given molecule is greater than 0.99 are used for the parameter learning for this given molecule;

• If the above criterion yields less than 5 molecules, we use 5 nearest molecules for the parameter learning. Here the 5-molecule nearest neighbor method is found to provide the best leave-one-out prediction.

Models for the total solvation free energy are given in Eq. (1). We utilize a unified protocol for the prediction of solvation free energy. First, for the polar solvation free energy ΔG^p, we select either the Poisson model with a given point charge force field or the polarizable Poisson model. Second, for the calculation of the non-polar solvation free energy, we note that in general, the functional group scoring approach can deliver better blind predictions than the nearest neighbor approach. Therefore, in our non-polar energy prediction step, the functional group scoring method is used whenever it works. Otherwise, we utilize the nearest neighbor approach. Our solvation prediction protocol is illustrated in a flowchart shown in Fig. 3.

We consider a total of 668 molecules, and the structures of these molecules are downloaded from the Pubchem project. The experimental solvation free energies of these molecules are collected from the literature. This dataset contains molecules from the Sampl blind solvation prediction projects, ranging from SAMPL0 to SAMPL4;^25,26,34,35 and the remaining molecules in our dataset are collected from the literature.^4,48,85 Coincidentally, there is a considerable overlap of our database with Mobley’s solvation database, which is available from http://mobleylab.org/resources.html. More precisely, 589 molecules in our database are already covered by Mobley’s. The detailed information of our dataset is provided in the supplementary material.

For both the standard Poisson model and the KSDFT based polarizable Poisson model, we consider four types of atomic radii, namely, Amber 6, Amber bondi, Amber mbondi2,⁵ and ZAP9⁶⁰ parametrizations. Additionally, for the standard Poisson model, three sets of charge assignments, namely, OpenEye-AM1-BCC v1 parameters,⁴⁰ Gasteiger,²⁴ and Mulliken,⁵ are tested. Our MIBPB solver^6,27,84 is utilized to solve the Poisson interface problem to obtain electrostatic solvation energies. The probe radius is set to 1.4 Å in all PB calculations.

For the KSDFT based polarizable Poisson model, the Poisson interface problem is coupled with the generalized KSDFT in a self-consistent approach. The charge density used by Poisson interface problem is calculated in ab initio approach by the generalized KSDFT, here generalization due to the additional solute-solvent reaction field energy is included in the effective KS potential. The Poisson interface problem is employed for calculating the reaction field energy. This approach is modified from our earlier differential geometry based polarizable PB model.⁹ We iteratively couple the SIESTA⁷⁵ and the MIBPB solver,^6,27,84 in which a sharp solvent-solute interface is employed. The SIESTA software with additional reaction field energy is used for charge density calculation, and our in-house software is used for reaction field energy calculation. A uniform mesh size of 0.25 Å is used for solving the Poisson equation in both the standard Poisson and Polarizable Poisson model. In the polarizable Poisson model, the computed change densities are mapped to the uniform mesh with the conservation of the total charge.

For the molecules containing iodine atoms, the current level of DFT method used in this work, including the Gaussian software, cannot handle this element appropriately. Therefore, for a uniform comparison, we ignore molecules containing iodine atoms. There is a similar situation for bromine atoms — there is no appropriate pseudo-potential for this atom in the SIESTA software.⁷⁵ Therefore, molecules involving bromine atoms are excluded in our KSDFT based polarizable PB calculations. Table III lists four molecules in SAMPLx that are not considered in some of our predictions.

For both the Poisson model and the polarizable Poisson model, the ESES software⁵⁰ is used for solute molecular surface generation. The MIBPB software^6,27,84 is used for solving the Poisson model. The polarizable Poisson model is solved by the self-consistent coupled MIBPB and SIESTA software.⁷⁵ 

In our non-polar solvation free energy model, both atomic surface areas and surface enclosed volume are required. In a recently developed Eulerian solvent excluded surface (ESES) package, a second-order accurate numerical scheme for surface area calculations and a third-order accurate numerical scheme for volume estimations have been developed.⁵⁰ A weighted Voronoi diagram algorithm is implemented to partition a molecular surface area into atomic surface areas. These schemes have been intensively validated by a large number of test examples. In this work, both atomic surface areas and molecular volume are calculated directly by using our ESES software package.⁵⁰ 

In this work, instead of using the classical non-polar model that includes the van der Waals interaction between the solvent and solute non-polar interaction⁸³ as shown in Eq. (7), we utilize the atomic surface area model given in Eq. (9) for the corresponding effects. In this part, we consider a set of numerical test to verify the validity of this treatment. Our numerical results indicate that the atomic surface area model provides results as good as those obtained by the van der Waals based non-polar model.

The van der Waals interaction Eq. (7) is modeled as a semi-continuous and semi-atomic potential. More specifically, we consider the 6-12 LJ potential for modeling the van der Waals interaction between the continuum solvent and atomistic solute molecules.⁸³ The probe radius of the solvent in the LJ potential is set to be 1.4 Å, which is the same as that used in all other calculations. Table IV lists the root mean square (RMS) errors with four types of atomic radii and four charge parametrization methods for the prediction of SAMPL0’s solvation free energies. From these results, it is seen that for a given force field parametrization, the atomic surface area approach provides slightly more accurate non-polar solvation free energy prediction. As shown in our earlier work,⁸³ the performance of the semi-continuous and semi-atomic LJ potential is sensitive to the choice of the probe radius, due to its nonlinear dependence on the radius, whereas the SES based atomic surface area approach is less sensitive to the probe radius. We conclude that in general the atomic surface area based non-polar model is at least as good as the LJ potential based one for modeling the non-polar solvation free energy. Therefore, the atomic surface area based non-polar model is employed in the rest of this paper.

We first examine the proposed nearest neighbor approach by the leave-one-out test of 668 molecules. In this examination, we select one molecule at a time and use all other molecules as the training set to predict the selected one’s solvation free energy. This process is applied systematically to all the molecules in the whole dataset of 668 molecules. Four different atomic radius sets are considered, together with four different charge force fields. Our results are illustrated in Fig. 4, and the associated values are listed in Table V. In general, all radius sets and charge force fields perform similarly well. The maximum RMS error is below 1.7 kcal/mol for all methods over all molecules. More specifically, Bondi and Mbondi radii offer the best overall results. For charge force field, AM1-BBC charges appear to provide the best predictions and their RMS errors are less than 1.5 kcal/mol for all the four radius sets. The optimal result obtained with AM1-BBC charges and ZAP9 radii has an RMS error of 1.33 kcal/mol. Figure 5 plots the optimal results on the leave-one-out test.

For a set of 643 molecules, which largely overlaps with present dataset, Mobley and Guthrie reported an RMS error of 1.51 kcal/mol.⁵⁷ Our results indicate that the present nearest neighbor approach can achieve highly accurate predictions for the solvation free energies for all the atomic radius and charge force fields.

In this section, we present results of the blind solvation free energy predictions based on the proposed protocol. All the SAMPL0-SAMPL4 challenges for solvation free energies are considered. We predict the solvation free energies with a leave-SAMPLx-out approach, in which the SAMPLx molecules and their experimental solvation free energies are regarded as unknown while information of other molecules is utilized to predict selected SAMPL test set, based on molecular formulas.

We first consider SAMPL0 test set, in which molecules are diverse. One molecule in this test set contains Br atom, for which our polarizable PB model does not an appropriate pseudo-potential. The structures of the SAMPL0 molecules can be found in the literature.⁶⁰ Tables S1-S4 of the supplementary material present the solvation free energies for the SAMPL0 test set obtained from both experiments and blind predictions using different charges and atomic radii parametrizations. The left chart of Fig. 6 shows the plot of RMS errors for each charge and atomic radius combination. It is clear that the change in atomic radii has a minor influence on the accuracy of predictions, while the change in charge force fields has a much more significant influence. In other words, the solvation free energy prediction for this test set is more sensitive to the charge parametrization than the atomic radii parametrization. In general, our KSDFT based polarizable Poisson model provides better predictions than those of other charge force fields. It is noted that the SIASTA based polarizable Poisson model with Amber Bondi radius parameters delivers the best solvation free energy prediction. This result is depicted in the right chart of Fig. 6. The associated RMS error (0.93 kcal/mol) appears to be better than that in the literature⁶⁰ (i.e., 1.71 kcal/mol for the full SAMPL0 test set of 17 molecules). The best prediction in the literature has an RMS error of 1.34 kcal/mol.⁴² Even though for the SIESTA based ab initio charge calculation, the molecule contains the Br atom is neglected, based on the general pattern that when any other force field used, the prediction for that omitted molecules is very well. We believe with proper pseudopotential for Br atom, our polarizable Poisson model can deliver an accurate prediction for this molecule also. Moreover, for this test set, all the other three charge parametrizations can provide a similar level of solvation free energy prediction.

We next consider the SAMPL1 challenge set for solvation predictions.³⁴ This set contains not only a largest number of 63 molecules, which is the largest among all the Sampl test sets, but also many molecules with extremely complicated structures. The detailed description of this set is given in the literature.³⁴ Most molecules in this set are druggable and very complex. The difficulty of this test set comes from two aspects: on one hand, the structures of molecules in this set are very complicated. On the other hand, the knowledge in the training set that can be used for predicting the solvation free energies of these molecules is rare or insufficient. In addition to the molecular complexity, the reported experimental solvation free energies also admit large uncertainty.³⁴ Most earlier computational predictions report RMS errors of 3–4 kcal/mol when some extremely complex molecules are excluded. Some best prediction for the whole set has an RMS error of 2.45 kcal/mol.⁴² The best performance was shown to give an RMS error of 2.4 kcal/mol on a subset of the SAMPL1 test set that contains only 56 molecules.⁵² Tables S5-S8 of the supplementary material present the solvation free energies for the SAMPL1 test set obtained from both experiments and blind predictions using different charges and atomic radii. Figure 7 plots the present blind predictions. One molecule (bromacil) that contains a Br atom is considered by all the charge models except for the polarizable PB method. It is noted that two of present approaches, namely, the AM1-BCC semi-empirical charges and the ab initio charge, provide relatively accurate predictions for this test set. When Amber Mondi2 atomic radii together with the SIESTA charge calculation is applied, the optimal solvation free energy prediction is achieved with RMS error 2.82 kcal/mol. When AM1-BCC charges and the Amber Bondi radii are used, the prediction for the whole SAMPL1 set without a molecule that contains Br atom has an RMS error of 3.06 kcal/mol, and the RMS error is 3.07 kcal/mol for the same test with all molecules. For this test set, the large prediction RMS error is mainly due to the extremely large error in predicting solvation free energies for some complex molecules, for which the RMS error can be as large as 15 kcal/mol. This unreasonable prediction is due to the fact that inappropriate molecule force field parametrization yields unreasonable electrostatic solvation free energies, which in turn leads to erroneous solvation free energy prediction. Similar to the SAMPL0 test set, the prediction is more sensitive to the charge parametrization. Nevertheless, the atomic radius parametrization also plays a critical role in the prediction accuracy for this test set.

SAMPL2 is another difficult test set with almost the same level of difficulty as the SAMPL1 test set.⁴⁴ Compared with SAMPL1 test set, it contains a few complex molecules, and most molecules in this test set are drug like ones as well. Contrary to molecules in SAMPL1 test set, this set has less uncertainty in the experimental solvation free energies. As shown in Tables S9-S12 of the supplementary material, the experimental solvation free energies of this test set distribute over a wide range. Using all-atom molecular dynamics simulations and multiple starting conformations for blind prediction, Klimovich and Mobley reported an RMS error of 2.82 kcal/mol over the whole set and 1.86 kcal/mol over all the molecules except several hydroxyl-rich compounds.⁴⁴ Some best prediction has an RMS error of 1.59 kcal/mol.⁴² In the present work, the molecule containing an I atom (5-iodouracil) is excluded in all calculations. Additionally, 5-bromouracil has a Br atom which is excluded in the polarizable PB model. Tables S9-S12 of the supplementary material list the details of the present predictions. The RMS errors from various radius and charge force fields are given Fig. 8. Apparently, this test set has a strong force field dependence as well. The RMS errors vary very much from one charge force field to another. However, the performance of these predictions has a weak radius dependence. The best prediction, obtained from a combination of Amber MBondi2 radius parameters and AM1-BCC charge force field, or Amber Bondi radius together with modified SIESTA charge, both have the RMS error of 1.90 kcal/mol when the molecule with a Br atom is excluded in the prediction. When all molecules are included, the combination of Amber MBondi2 radius and AM1-BCC charge parametrization gives the optimal prediction of RMS error 1.96 kcal/mol as depicted in the right chart of Fig. 8. Compared to the prediction of SAMPL1 test set, these predictions are more accurate, due to two reasons. First, the molecules in this set are slightly simpler and experimental uncertainly is less severe for the deterministic prediction. Second, in our parametrization of the non-polar solvation free energy prediction, we have more similar molecules in training set, which enables us to obtain better non-polar solvation free energy parameters in the nearest neighbor based approach. In contrast, in the SAMPL1 prediction, when the nearest neighbor approach is applied, the nearest molecules selected from the training set differ much from SAMPL1 molecules.

SAMPL3 test set is a relatively easy one with 36 solute molecules.²⁵ Its molecular structures are less versatile than the earlier test sets. Additionally, in this test set, there is no molecule that involves Br or I atom. One of difficulties in the prediction of this set is the lack of similar molecules in our database when all the SAMPL3 molecules are left out. The best prediction in the literature offers an RMS error of 1.29 kcal/mol.⁴² Our blind predictions and experimental results for SAMPL3 test set are given in Tables S13-S16 of the supplementary material. The RMS errors of our blind predictions for the SAMPL3 test set are plotted in Fig. 9. In this case, the accuracy of predictions also depends strongly on charge force fields and weakly on radius parameters. In general, Gasteiger charges are superb for this test set and their RMS errors are always smaller than 1 kcal/mol. Our best prediction, obtained from the combination of ZAP9 radius parameters and Gasteiger charge force field, has a small RMS error of 0.78 kcal/mol and is depicted in the right chart of Fig. 9. In general Gasteiger charge performs very poor in the solvation free energy prediction for previous tests on complex molecules, while we see that this charge parametrization method is superior for chlorinated hydrocarbon molecules in the present test set. Unfortunately, there is no uniformly optimal parametrization for all the molecules. This fact motivates us to seek a solvation free energy prediction that does not heavily depend on the force fields parametrization.

SAMPL4 test set is studied by using the proposed methods as well. A key feature of this test set is that its molecules are diversified. As shown in Tables S17-S20 of the supplementary material, this test set also involves a wide range of solvation free energies. However, the structures of these molecules are not as complex as those in both SAMPL1 and SAMPL2 test sets, which indicates a slightly easier task for solvation free energy prediction. This test set was studied by a number of researchers in the literature and the best prediction in the literature has the RMS error of 1.2 kcal/mol.⁵⁹ Figure 10 illustrates the RMS errors of our predictions for a total of 16 charge and radius combinations. Compared with those in the literature, all of our predictions are of high quality. Our best result has the RMS error of 1.03 kcal/mol and the corresponding results are given in the right chart of Fig. 10.

Table VI lists the root mean square errors (RMSEs), mean errors (MEs), Pearson correlations (Rs), and P-values (Ps) of our blind predictions of solvation free energies for all five test sets using a total of 16 charge and radius implementations. Some comments are in order. First, all predictions are quite sensitive to charge force fields and relatively less dependent on radius parameter selection. For SAMPL0 test set, although predictions with the AM1-BCC charge set have similar correlations as those obtained with other charge force fields, their MEs are significantly larger. However, AM1-BCC charge force field does a relatively better job in SAMPL1 and SAMPL2. Additionally, it is difficult to identify a clear winner over all the test sets — some approaches perform better in one or two test sets, but do not do well in the rest test sets. This phenomenon highlights the difficulty of designing optimal models for solvation analysis. Moreover, the KSDFT based polarizable Poisson model is more often to provide better predictions over most charge force fields and radius parameters. This indicates a potential to develop better physical models by improving the quantum charge density calculation. Further, the Pearson correlation becomes an interesting measure in SAMPL3 test set. Results obtained with AM1-BCC and Mulliken charge sets have particularly low correlations, although their RMSEs are quite good. Results for SAMPL1 and SAMPL4 test sets are quite consistent in terms of their RMSE, ME, and correlations. Furthermore, our P-value analysis indicates that most results very reliable, except for cases when the AM1-BCC and Mulliken charge forces fields are applied to SAMPL3. However, all P-values are smaller than the standard cutoff, 0.05. Finally, we point out that the blind prediction results presented in the present work are the state-of-the-art compared to those in the literature.^{21,25,26,34,35,42,55,56,58–60,68,69}

In this work, we propose a protocol to parametrize implicit solvent models for the blind prediction of solvation free energies. The proposed protocol parametrizes non-polar solvation model by scoring and ranking algorithms while utilizes either the Poisson model or a density functional theory (DFT) based polarizable Poisson model for polar solvation free energy calculations. For the parametrization of non-polar solvation models, we first utilized the assumption that molecules with the same functional group admit the same parametrization of the non-polar solvation energy functional. For complex poly-functional-group molecules, we develop a scoring procedure to determine the optimal relative weight of each functional group. For extremely complex molecules that fail the functional group scoring method, we further develop a molecule ranking algorithm to select an optimal set of nearest neighbor molecules for parameter training. We construct atomic features for the molecule ranking. Finally, we systematically integrate the above mentioned models and algorithms into a robust protocol for blind solvation free energy prediction.

In the present work, we considered an experimental database of 668 solvation molecules, see Table S21 in the supplementary material, the largest database ever constructed for solvation, to validate our approach. Among them, SAMPL0–SAMPL4 test sets are paid special attention. For the Poisson model or DFT based polarizable Poisson model, four sets of atomic radius parameters (i.e., Amber 6, Amber bondi, Amber mbondi, and ZAP9 radii) are combined with four sets of charge force fields (i.e., AM1-BCC, Mulliken, Gasteiger, and SIESTA DFT) to arrive at a total of 16 different implementations. The resulting polar solvation free energies are utilized in our parametrization for blind predictions. We first carry out the leave-one-out validation of the whole database. The AM1-BCC charge force field delivers a low root mean square (RMS) error of 1.33 kcal/mol, which is the lowest for such a large test database, to our best knowledge. We further conduct a series of leave-SAMPLx-out blind tests. On average, the AM1-BCC parametrization in the Poisson model and DFT based polarizable Poisson model performs better than other charge force fields, especially for predicting the solvation free energies of the complex molecules. The optimal RMS errors for SAMPL0-SAMPL4 are, respectively, 0.93, 2.82, 1.90, 0.78, and 1.03 kcal/mol.

From the solvation free energy predictions, particularly on SAMPL1 and SAMPL2 test sets, we conclude that atomic charge parametrization is extremely important for the present physical models, namely, the Poisson model or the KSDFT based polarizable Poisson model. Without an appropriate charge parametrization, the prediction errors can be amplified for molecules with complex structures. In general, for four charge parametrization methods, both the semi-empirical AM1-BCC charge and the ab initio charge calculation from the generalized KSDFT can provide relatively reliable charge assignments. For this reason, a solvation free energy prediction method that does not heavily depend on the molecule parametrization is under our consideration. Essential idea is that if one does not partition the solvation free energy into two isolated parts, the prediction errors in the polar solvation free energy will not propagate to the non-polar solvation free energy prediction. Alternative, our differential geometry based solvation model that dynamically couples polar and non-polar models^7,89 might also provide a less sensitive approach.

This work can be improved in a number of ways. First, the classification of the database into functional groups is not unique. Future study will explore optimal molecules partition. Additionally, the selection and computation of atomic features need to be further investigated. It is possible to construct an optimal set of atomic features for solvation analysis and prediction. Further, in the current work, molecular ranking for nearest neighbor searching is not optimal yet. More sophisticated machine learning and/or deep learning algorithms can be developed for this purpose. Finally, a more versatile DFT solvers can be utilized to further improve our in-house polarizable Poisson model. These issues are under our consideration.

See supplementary material for specific predictions of SAMPL0-SAMPL4 test sets for all force field parametrizations and molecular name, PubChem ID, and experimental solvation free energies of all 668 molecules. Additionally, more detailed datasets with coordinates, radii, and charges for 668 molecules can be downloaded from website http://weilab.math.msu.edu/MIBPB/.

This work was supported in part by NSF Grant No. IIS-1302285 and MSU Center for Mathematical Molecular Biosciences Initiative. B.W. thanks Dr. Pengfei Li, Professor Ray Luo, Professor David Mobeley, and Professor Wei Yang for helpful discussions. G.-W.W. acknowledges Nathan Baker, Michael Gilson, and Weitao Yang for useful discussions.