Implicit Solvent Models and Their Applications in Biophysics

1. Introduction

Solvation phenomena profoundly influence the structure, dynamics, and function of biomolecules in aqueous and non-aqueous environments [1]. In biological systems, water and other solvents mediate essential processes such as protein folding, ligand binding, molecular recognition, and catalysis [2,3]. Accurately modeling solvation is important in computational biophysics and chemistry with implications ranging from fundamental molecular science to practical drug discovery and (bio)materials design [4]. Traditionally, explicit solvent models, where each solvent molecule is treated as a discrete particle, have been the gold standard for capturing solvation effects in molecular simulations. While explicit approaches provide detailed insights into solvent structure and dynamics, their high computational cost and the need for extensive sampling have motivated the development of alternative strategies [5]. Among these, implicit solvent models have emerged as crucial tools, offering a balance between computational efficiency and physical realism by replacing discrete solvent molecules with a continuum representation [6].

The conceptual foundations of implicit solvent models can be traced to early dielectric theories of solvation. The seminal work of Onsager and Debye in the early 20th century established the treatment of solvents as dielectric continua, enabling the estimation of solvation energies based on bulk properties such as dielectric constant and molecular polarizability [7]. These early models laid the groundwork for the development of more sophisticated theoretical frameworks [8]. With the advent of computational chemistry, the need for practical and generalizable solvation models led to the development of continuum electrostatic approaches[9]. The Poisson-Boltzmann (PB) equation provided a rigorous description of electrostatic interactions between solutes and a surrounding dielectric medium, representing spatial variations in dielectric properties and ionic strength [9,10]. Furthermore, the Generalized Born (GB) model introduced efficient pairwise approximations to the PB formalism, enabling rapid estimation of electrostatic solvation energies for large biomolecular systems [9,11]. In parallel, advancements in quantum chemistry drove the creation of implicit solvation models such as the Polarizable Continuum Model (PCM) and the Conductor-like Screening Model (COSMO), which facilitated the inclusion of solvation effects in electronic structure calculations [12,13]. More recent models, such as the SMx family and SMD, integrate both electrostatic and non-electrostatic (cavitation, dispersion, and repulsion) contributions to provide highly accurate predictions of solvation free energies across a wide range of solvents and solutes [12,14,15].

At the center of implicit solvent models lies the partitioning of the solvation free energy into physically meaningful components [9]. Typically, this includes a polar (electrostatic) term, accounting for the interaction of the solute’s charge distribution with the dielectric environment, and a nonpolar term, describing contributions from cavity formation, solvent-accessible surface area, and van der Waals interactions [9,16,17,18]. These terms have been refined to incorporate additional physical effects, such as hydrogen bonding corrections and solvent-excluded volume [17]. Implicit solvent models are grounded in continuum theories, where the solute is embedded in a dielectric medium characterized by macroscopic properties. The electrostatic component is frequently computed by solving the Poisson-Boltzmann equation or its linearized form, or by employing the Generalized Born approximation. The nonpolar component is often related to the solvent-accessible surface area (SASA) or volume (SAV), with empirical parameters derived from experimental data or explicit solvent simulations [18,19,20,21]. Modern models further decompose the nonpolar term to distinguish between repulsive (cavity) and attractive (dispersion) interactions [18,22,23]. The primary advantage of implicit solvent models lies in their computational efficiency [9,18,24]. By removing the need to simulate thousands of explicit solvent molecules, these models enable rapid exploration of biomolecular conformations, facilitate enhanced sampling, and make feasible the simulation of large or complex systems that would be otherwise impossible. Thus, implicit solvent models have been widely adopted in biomolecular simulations, structure-based drug design, protein-protein and protein-ligand binding studies, and nucleic acid research [9]. However, the continuum approximation introduces inherent limitations. The absence of explicit solvent structure can lead to inaccuracies in capturing specific solvent-mediated interactions, such as water bridges, hydrogen bonds, and ion effects. Implicit models may also struggle to represent the entropic contributions of solvent molecules and the heterogeneous nature of biological environments [5,25]. Parameterization remains a challenge, as the accuracy of these models depends strongly on the choice of atomic radii, dielectric constants, and empirical coefficients [26]. Despite these challenges, ongoing refinement and the integration of machine learning (ML) techniques continue to enhance the reliability and predictive power of implicit solvent approaches [27].

The versatility of implicit solvent models has catalyzed their adoption in a wide array of biophysical applications [9]. In protein-ligand binding, implicit solvation methods are routinely employed to estimate binding free energies, rank inhibitor potency, and guide lead optimization [9]. For intrinsically disordered proteins (IDPs) whose lack of stable tertiary structure challenges traditional modeling, implicit solvents enable the efficient exploration of vast conformational landscapes [28,29,30] and facilitate comparison with experimental data such as FRET or SAXS [31,32]. In nucleic acid simulations, hybrid explicit and implicit solvent approaches have yielded accurate predictions of DNA and RNA structure and dynamics with significant reductions in computational cost. Moreover, advances in quantum implicit solvent models allow for the incorporation of solvation effects in ab initio calculations of reaction mechanisms, spectroscopic properties, and electronic structure predictions [12]. This has implications for the study of enzyme catalysis, photochemistry, and (bio)materials design [12,33].

Here, we provide an overview of implicit solvent models, detailing their theoretical foundations, practical implementations, and applications in biophysics and computational chemistry. We describe classical and quantum formulations, benchmark studies, and recent developments in model parameterization and hybrid approaches. Special attention is given to recent ML-augmented implicit solvent models and their impact on accuracy and transferability. We also discuss current challenges and highlight future directions, including the integration of quantum computing and multi-scale modeling. By elucidating both the strengths and limitations of implicit solvent models, we aim to provide a critical resource for the selection and application of appropriate solvation models in diverse biophysical and (bio)chemical investigations.

2. The Beginning of an Era

Onsager model was introduced by Lars Onsager [34] to refine Lorentz local field approximation in a system of dipoles [35]. Debye’s dielectric constant,$\epsilon$, equation is given in the following:

$${\displaystyle \frac{\epsilon -1}{\epsilon +2}}={\displaystyle \frac{4\pi }{3}}\sum N(\alpha +{\displaystyle \frac{{\mu }^2}{3kT}})$$

(1)

where $N$ is the concentration of molecule, $\alpha$ is the polarizability of a molecule, and $\mu$ is molecule’s permanent dielectric moment [34]. Debye’s equation failed in predicting pure polar liquids when compared with experiments. In Kirkwood’s work,³ it was observed that for dilute solutions of polar substances in non-polar liquids, the dipole moments calculated using the theoretical formula closely match those measured in the vapor phase, though there is a systematic trend for the moment (p) to decrease as the dielectric constant of the environment increases. However, for pure polar liquids, the calculated dipole moments are significantly smaller, and this discrepancy becomes more pronounced as the dielectric constant of the liquid rises.

The solvation free energy, $∆G_{solv}$ although its definition and content have changed over time [9,16,17,18,24,36], plays a vital role in for thermodynamical, physical description of a solution [9,24], but in order to explain the free energy of solvation in a comprehensible way, these differences will be discussed in detail. Hereby, we first introduce nonpolar free energy of solvation and then the other approaches that may or may not incorporate with nonpolar free energy component. Nonetheless, $∆G_{solv}$ is composed of two main components [5,16,36]: polar and apolar components that are denoted by different appellations [5,17,36,37]. Electrostatic free energy of solvation, $∆G_{ele},$ is used to refer to electrostatic interactions between solute and solvent [9,17,24,36,37], whereas $∆G_{np}$ denotes the nonpolar interactions:

$$∆G_{solv}=∆G_{np}+∆G_{ele}$$

(2)

In parallel, $∆G_{solv}$can be separated into three main components [9,24,38,39]. Firstly, $∆G_{cav}$ is used to define the excluded volume, also cavity for the solvent around the solute [9,24,39,40]; $∆G_{ele}$ is used to refer to electrostatic interactions between solute and solvent [9,24]; and finally $∆G_{vdW}$ describes the van der Waals interactions [9,24]:

$$∆G_{solv}=∆G_{cav}+∆G_{ele}+∆G_{vdW}$$

(3)

Tan et. al. stated that the nonpolar component of the solvation free energy $∆G_{np}$ consisted of two components of highly different interactions, which are repulsive free energy and attractive free energy and therefore must be modeled separately [18]. Gonçalves et. al. introduced a short range free energy instead of crude solvation free energy of the van der Waals interactions $∆G_{vdW}$, which incorporates small-scale attractive and repulsive forces between solute and solvent molecules [24]. Before, it was thought that for transferring a molecule from the gas phase to the solvent, the cavity can be treated as spherical as referred in equation 4 [41]:

$$∆G_{solv}=4\pi r^2\sigma -ϵ,$$

(4)

where $r$ is the cavity radius, $\sigma$ is the solvent surface tension, and $ϵ$ is the solute-solvent interaction energy [41]. However as mentioned above, the solvation free energy has polar and nonpolar components [17,20], yet these components and later their sub-components would only be relevant in further developments of implicit solvent models [9,16,17,18,20,42]. Next, accessible surface area (ASA) was defined[43] where ASA refers to the sum of the area around the atom of a molecule - that solvent molecules can interact with - without penetrating any other atoms of the molecule. For calculating ASA, different methods have been introduced [44,45,46]. However, ASA calculations required that buried atoms to be considered as they could not contribute to the free energy of solvation [46,47]. ASA for an atom $A_i$ and its atomic solvation parameter $∆{\sigma }_i$ of the free energy of transfer is (aqueous solution) [47,48]:

$$∆G_i=∆{\sigma }_i{A}_i$$

(5)

where the sum of free energy transfer for a whole molecule of interest $∆G_R$, such as a residue of a protein can be described as:

$$∆G_R={\sum }_i∆{\sigma }_i{A}_i$$

(6)

However, we have to mention that polar interactions are omitted and the calculation solely depends on ASA. Moreover, atomic solvation parameters vary for each atom subtypes. Thus, ASA parametrizations were then conducted [9,47,48]. Later, by Stil et al., it was proposed that nonpolar solvation, cavity and van der Waals interactions can be defined as a function of ASA [49]:

$$G_{cav}+G_{vdW}=\sum {\sigma }_{k}SA{S}_k$$

(7)

where $SA{S}_k$ is the total solvent-accessible surface area of atoms of type $k$ and ${\sigma }_k$ is an empirical atomic solvation parameter. $∆G_{np}$ is a key component of $∆G_{solv}$, thus, nonpolar solvation free energy can be defined in terms of SASA [18]:

$$∆G_{np}=\gamma SAS+c$$

(8)

where $\gamma$ represents, which is surface tension coefficient, nonpolar energy representation per unit surface area and $c$ is a constant where it defines nonpolar solvation energy. However, it has been addressed that, for small hydrophobic solutes, which have protein-like functional groups, solvation free energy approximation is better when the volume of the solute instead of the surface area was used [50,51]. Moreover, it was found that - on the fine grained energy scale necessary to predict the high-resolution structure of proteins and protein-ligand complexes - the correlation between $∆G_{vdW}$ and the SASA of the solute is poor [16,17,18]. Therefore, modeling $∆G_{vdW}$ and $∆G_{cav}$ independently is better for accurate representation. Levy et. al. proposed also the backbone of the AGBNP model [16,17]. Since $∆G_{vdW}$ and $∆G_{cav}$ are defined separately in the literature; $∆G_{cav}$ can be described as [22,52,53,54,55,56,57,58]:

$$∆G_{cav}= \gamma SAS$$

(9)

where γ represents surface tension coefficient.

Levy et. al. [17] introduced more specific parameters and definitions for this model:

$$∆G_{cav}={\sum }_i{\gamma }_i{A}_i$$

(10)

where $\gamma$ defines surface tension parameters of the corresponding atom, $A_i$ denotes the van der Waals surface area of corresponding atom and the summation is over solute atoms. However, in this model, instead of the free energy of cavity formation based on SASA, another approach based on van der Waals surface area and surface tension parameters was used [16,17,20]. They also described the solute–solvent van der Waals free energy, $∆G_{vdW}$, by an expression, which was obtained as an integral of the van der Waals solute–solvent interactions over the solvent volume modeled, as a uniform continuum [16,17,20]:

$$∆G_{vdW}={\sum }_i{\alpha }_i{\displaystyle \frac{a_i}{{\left(B_i+R_w\right)}^3}}$$

(11)

where $B_i$ is the Born radius of $i$’th atom, $R_w$ is the radius of the water molecule, ${\alpha }_i$, which is a scaling factor and given by [16,17,20]:

$${\alpha }_i=-{\displaystyle \frac{16}{3}}\pi {\rho }_w{\epsilon }_w{\sigma }_{iw}^6$$

(12)

${\rho }_w$ is the density of the water molecule, ${\sigma }_{iw}$ and ${\epsilon }_w$ are the parameters of the OPLS force field [20,21].

As previously mentioned, the original AGBNP model relied solely on the van der Waals surface, calculated using the Poincaré formula [59], to estimate the solute’s surface area and volume [16,17]. This approach could lead to inaccurate approximations of Born radii and, consequently, did not adequately represent the exposed surface area in the model. To address these shortcomings, Levy et al. developed the AGBNP2 model, introducing a new term: the hydrogen bonding correction. The total solvation free energy in AGBNP2 is given by [16]:

$$G_{solv}=∆G_{cav}+∆G_{ele}+∆G_{vdW}+∆G_{hb}$$

(13)

where $∆G_{hb}$ equals to:

$$∆G_{hb}={\sum }_s{h}_{s}S(w_s;w_a,w_b)$$

(14)

where $h_s$ is the maximum correction energy which depends on the type of solute-water hydrogen bond and $S(w_s;w_a,w_b)$ is a polynomial switching function, which can be found in detail in Ref. [16].

Later, it was proposed that nonpolar solvation ($∆G_{np})$ can be defined with repulsive free energy $(∆G_{rep})$ and attractive free energy $(∆G_{att})$ sub-components [18]. The nonpolar repulsive free energy ($∆G_{rep})$ can be modeled as [18]:

$$∆G_{rep}= \gamma SAS+c$$

(15)

However, it has been observed that solvation free energy approximation depends on the usage of solvent-accessible volume or area. The applicability of SASA for describing hydrophobic solvation is limited to larger length scales; for surfaces with significant roughness or high curvature, solvent-accessible volume (SAV) serves as a more appropriate metric. It has been shown that for small cavities, repulsive free energy $(∆G_{rep})$ correlates with the volume of the cavity, whereas for large cavities, it correlates with the cavity surface [18,19]:

$$∆G_{rep}= pSAV+c$$

(16)

where SAV refers to solvent-accessible volume and $p$ is a solvent pressure parameter.

On the other hand, Chandler et.al. [18,22,23] showed that the attractive free energy $(∆G_{att})$ can be described as [18,22,23]:

$$∆G_{att}={<U}_{att}^{uv}>$$

(17)

where ${<U}_{att}^{uv}>$ is the van der Waals attractive interaction potential energy between solute $u$ and solvent $v$ and described as the solvent-occupied volume integration [18]:

$${<U}_{att}^{uv}>= {\sum }_{a=1}^{N_s}\int {\rho }_{aw}\left({\mathbf{r}}_{\mathrm{a}\mathrm{w}}\right){\mathrm{V}}_{\mathrm{a}\mathrm{t}\mathrm{t}}\left({\mathbf{r}}_{\mathrm{a}\mathrm{w}}\right)d{\mathbf{r}}_{\mathrm{a}\mathrm{w}}$$

(18)

where $N_s$ is the sum over all solute atoms and ${\rho }_{aw}\left({\mathbf{r}}_{\mathrm{a}\mathrm{w}}\right)$is a solvent distribution function around the solute atom $a$ at a given solute-solvent distance, ${\mathbf{r}}_{\mathrm{a}\mathrm{w}}$; ${\mathrm{V}}_{\mathrm{a}\mathrm{t}\mathrm{t}}\left({\mathbf{r}}_{\mathrm{a}\mathrm{w}}\right)$ is the attractive van der Waals potential in a decomposition scheme. However, due to the impossibility of knowing the value for ${\rho }_{aw}\left({\mathbf{r}}_{\mathrm{a}\mathrm{w}}\right)$ without equilibrium simulations in explicit solvent, it must be approximated [18].

3. Poisson-Boltzmann Equation

As mentioned above, solvation free energy is acquainted with electrostatic properties, $∆G_{ele}$, that has to be calculated and it was proposed that it can be calculated by solving the Poisson-Boltzmann equation with the continuum dielectric model [18,49,60]. The Poisson-Boltzmann equation is a widely used equation to describe the behavior of charged molecules in ionic solutions [10,61,62,63,64]. This equation allows to model electrostatic interactions by combining the principles of Poisson equation and Boltzmann distribution [9,10,61,62,63,64]. First attempts to solve it involved basic geometric models such as spheres for proteins [65] and cylinders for DNA [66]. Later, more accurate representations of the molecular surface were developed [67], such as the Solvent Excluded Surface [68], the Minimal Molecular Surface [69], and, more recently, the promising Blobby Surfaces, which utilize a multi-level Gaussian density map [70]. Nevertheless, the Poisson equation relates the scalar electric potential (Φ) to the charge density (ρ) [9]. The Poisson-Boltzmann equation can be expressed as follows [9,71]:

$$\overrightarrow{\nabla }\left[\epsilon \left(\mathit{r}\right)\overrightarrow{\nabla }\Phi \left(\mathit{r}\right)\right]=-{\rho }_B\left(\mathit{r}\right)-{\sum }_i{c}_i^b{z}_{i}q{e}^{-z_{i}q\Phi \left(\mathit{r}\right)*\beta }$$

(19)

where $\epsilon$ is the local dielectric constant, $c_i$ is the bulk concentration of ion i, $z_i$ is the valence, $q$ is the unit charge, $\beta$, which equals to ${\displaystyle \frac{1}{k_{B}T}}$. In such cases, when the exponential term is sufficiently small, it can be approximated by linearization. This leads to the formulation of the Linear Poisson-Boltzmann equation:

$$\overrightarrow{\nabla }\left[\epsilon \left(\mathit{r}\right)\overrightarrow{\nabla }\Phi \left(\mathit{r}\right)\right]=-{\rho }_B\left(\mathit{r}\right)-{\sum }_i{c}_i^b{z}_i^2{q}^2{\displaystyle \frac{q^2\Phi \left(\mathit{r}\right)}{1/\beta }}$$

(20)

Notably, Poisson-Boltzmann Equation and its discretization methods are explained in detail in the literature [72,73]. These are used to determine electrostatic free energy and can be shown as a function of the electrostatic part of the solvation free energy:[74]

$$∆G_{ele}={\displaystyle \frac{1}{2}}\int {\Phi }_{reac}\left(\mathit{r}\right)\rho \left(\mathit{r}\right)dV$$

(21)

where ${\Phi }_{reac}$ is the reaction field which is the difference between the potentials which describe the solvent ${(\Phi }_{sol})$and the vacuum${(\Phi }_{vac})$.

Electrostatic free energy can also be decomposed into individual components. For example, Honig et al. have provided a detailed description of this decomposition in their studies, which will be briefly summarized here [36,60,75]:

$$∆G_{ele}=∆G_{ele}^{fixed}+∆G_{ele}^{ions}+∆G_{ele}^{sol}$$

(22)

where $∆G_{ele}^{fixed}$ is the fixed charge of solute, which equals to:

$$\mathsf{\Delta }{G}_{ele}^{fixed}={\displaystyle \frac{1}{2}}{\int }_V{\rho }^{\mathit{fixed}}(\mathit{r})\phi (\mathit{r})dv$$

$$\mathsf{\Delta }{G}_{ele}^{ions}=k_B\mathrm{T}{\int }_V{\sum }_i[{\mathrm{c}}_{\mathrm{i}}^{\infty }-{\mathrm{c}}_{\mathrm{i}}(\mathit{r})]dv$$

$$\mathsf{\Delta }{G}_{ele}^{sol}=-{\displaystyle \frac{1}{2}}{\int }_V{\rho }^{\mathit{sol}}(\mathit{r})\phi (\mathit{r})dv$$

Later, Nguyen et. al. showed that electrostatic solvation free energy $(∆G_{ele})$ can depend on Poisson-Boltzmann:

$$∆G_{ele}={\displaystyle \frac{1}{2}}{\sum }_{i=1}^{N_m}{q}_i(\mathsf{\Phi }({\mathbf{r}}_{\mathrm{i}})-{\mathsf{\Phi }}_0({\mathbf{r}}_{\mathrm{i}}))$$

(23)

Here, Φ(r_i) is the electrostatic potential at the position vector of the ith particle, Φ0​ is the solution of the Poisson-Boltzmann equation without accounting for the solvent–solute interface, and qiis is the partial charge of the ith atom located at position r_i.

4. Born Equation

Born equation describes the transfer free energy of a spherical ion from the gas phase to a dielectric medium, where both the potential and solvation free energy are determined based on [11,36,74,76]:

$$\mathsf{\Delta }{G}_{ele}=-{\displaystyle \frac{q^2}{2a}}(1-{\displaystyle \frac{1}{{\epsilon }_w}})$$

(24)

where $a$ is the ion radius, $q$ is its charge and ${\epsilon }_w$ is the continuum dielectric, which represents water. It can also be generalized as [9,20,77]:

$$\mathsf{\Delta }{G}_{ele}=-{\displaystyle \frac{q_i^2}{R_i}}({\displaystyle \frac{1}{{\epsilon }_{in}}}-{\displaystyle \frac{1}{{\epsilon }_{out}}})$$

(25)

where $R_i$ represents the Born radii.

The electrostatic solvation free energy, depending on the green function via the generalized Born model, exists in different forms in the literature [36,49]. The electrostatic free energy part of the solvation free energy can be formulated as [11,36]:

$$∆G_{ele}=-{\displaystyle \frac{1}{2}}(1-{\displaystyle \frac{1}{{\epsilon }_w}}){\sum }_{i,j }^N{\displaystyle \frac{q_i{q}_j}{f_{GB}}}$$

(26)

where $f_{GB}$is a function that relates the effective Born radii and has been defined in the literature [9,11,36,49,74,78]. Notably one of the reliable Generalized Born formula was proposed by Still et. al. [77]. Sigalov et. al. [76] expanded the formulation beyond the canonical form of the green function:

$$∆G_{ele}=-{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{{ϵ}_{in}}}-{\displaystyle \frac{1}{{ϵ}_{out}}}\right){\displaystyle \frac{1}{1+\beta \alpha }}{\sum }_{ij}{q}_i{q}_j ({\displaystyle \frac{1}{f^{GB}}}+{\displaystyle \frac{\alpha \beta }{A}})$$

(27)

Here, $\beta$ is the ratio of the dielectric constants, A is the electrostatic size of the considered molecule, and $\alpha$ is a constant parameter.

Lazardis et. al. proposed another model (EEF1), which is a less heuristic approach and combines an excluded volume approach with a modified version of the polar hydrogen energy function [79]. Since the Solvation Free Energy$, ∆G_{solv},$ of a given macromolecular conformation can be written as an integral over the space around it:

$$∆G_{solv}={\sum }_i∆G_{solv}^i$$

(28)

Moreover, it can be represented as an integral over the space around it:

$$∆G_{solv}={\sum }_i\int f\left(\mathit{r}\right)d\mathit{r}$$

(29)

$f\left(\mathit{r}\right)$ is the solvation free energy density at point $\mathit{r}$. This can be further improved by treating the solute-solvent interaction energy and using the correlation function:

$$∆G_{solv}^i=∆G_{ref}^i-{\sum }_j{\int }_{Vj}{f}_i\left(\mathbf{r}\right)\mathrm{d}\mathbf{r}$$

(30)

This can be decomposed into a sum of pairwise interaction:

$$∆G_{solv}={\sum }_i∆G_{solv}^i={\sum }_i∆G_{ref}^i-{\sum }_i{\sum }_{i\ne j }{f}_i\left(r_{ij}\right)V_j$$

(31)

Here, $∆G_{ref}^i$ refers to the reference Solvation Free Energy in which $i$ is essentially fully solvent-exposed group [79]:

$${\sum }_i∆G_{ref}^i{\sum }_{j>i }{(f}_i\left(r_{ij}\right)V_j+f_j\left(r_{ij}\right)V_i)$$

(32)

Thus,

$${\sum }_i∆G_{ref}^i-{\sum }_{j>i}\{{\displaystyle \frac{2∆G_{free}^i}{4\pi \sqrt{\pi }{\lambda }_i{r}_{ij}^2}}{e}^{-x_{ij}^2}{V}_j+{\sum }_{j>i}\{{\displaystyle \frac{2∆G_{free}^i}{4\pi \sqrt{\pi }{\lambda }_i{r}_{ij}^2}}{e}^{-x_{ji}^2}{V}_j \}$$

(33)

where $∆G_{free}^i$ refers to the solvation free energy of the isolated group and ${\lambda }_i$ is a correlation length [79].

5. The DelPhI Model

Developed in the early 1990s and acquired many updates since then, DelPhI is designed to solve the Poisson-Boltzmann equation for biomolecular electrostatics and it is still widely used [67,80,81,82,83]. Due to its continuous expansion, we will introduce its fundamentals. Here, we introduce a version of the Poisson-Boltzmann Equation:

$$\nabla \left[\epsilon \left(\mathit{r}\right)\nabla \Phi \left(\mathit{r}\right)\right]-K{\left(x\right)}^2\sinh \left(\Phi \left(x\right)\right)=-4\pi \rho (x)$$

(34)

where $\Phi \left(x\right)$ is the electrostatic potential and is determined by $\epsilon \left(x\right)$, the spatial dielectric function, $K{\left(x\right)}^2$,which is a modified Debye-Hückel parameter, and reason lies behind that is the dielectric discontinuity resulting from the difference in the electronic polarizability between the macromolecule and the surrounding solvent, which significantly affects electrostatic interactions, and lastly, $\rho (x)$ is the charge distribution function[80]. Linearization is carried out by approximating the hyperbolic sine function using its argument, effectively replacing the non-linear $\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}$ term with a linear expression [80].

They also considered cubic lattices, in which each side contains $L$ grid points, therefore total of $N=L^3$ points. Thus, the final reduction to finite difference form yields:

$${\Phi }_0=\left[{\displaystyle \frac{\left({\sum }_{i=1}^6{\epsilon }_i{\Phi }_i\right)+4\pi q_0/h }{\left({\sum }_{i=1}^6{\epsilon }_i\right)+{\left(K_{0}h\right)}^2}}\right]$$

(35)

where ${\Phi }_0$ denotes the potential at a particular grid point, ${\Phi }_i$ the potential at the six nearest neighbors, ${\epsilon }_i$the dielectric constant at the midpoint between ${\Phi }_0$, and ${\Phi }_0$, the charge assigned to the grid point is given by $q_0$, and $K_0$ is nonzero if the grid point is in salt. The parameter $h$ denotes the grid spacing [80].

Also, it can be written as a matrix:

$$\Phi =\mathit{T}\Phi +\mathit{Q}$$

(36)

where $\mathit{T}$ denotes a matrix, $\Phi$ and $\mathit{Q}$ are vectors. The elements of the matrix are:

$$T_{ij}=0 if j\ne g(x\pm 1,y,z),g(x,y\pm 1,z),g(x,y,z\pm 1)$$

(37)

And the elements of $\mathit{Q}$ are:

$$Q_i=[4\pi q(x,y,z)/h]*d(x,y,z)$$

(38)

However, boundary conditions are crucial and are typically fixed by using Coulomb's law or Debye-Hückel theory. At a sufficient distance from the molecular surface, the potential becomes insensitive to the specific molecular shape, thus making the choice of box fill ratio important[ 80].

Additionally, Jacobian relaxation method is prone to slow convergence as the grid increases thus making large-scale computations costly both in terms of time and memory. To address these issues Gauss-Siedel relaxation can be employed, which uses newly updated potential values during the iteration process. Unlike Jacobian, where the potential at a current iteration is calculated solely based on the potential of the previous iteration, Gauss-Seidel incorporates recent potential values into the ongoing iteration. However, the implementation order of grid points can vary, thus the spectral radius and intrinsic rate of convergence of the method remain independent of this mapping. Due to nature of the matrix which is consistently ordered, for the Poisson–Boltzmann problem, Gauss–Seidel converges approximately twice as fast as the Jacobian method, requiring only half as many iterations [80] which leads to the spectral radius of Gauss-Seidel being the square of that for Jacobian:

$${\lambda }_{N}Siedel={\lambda }_N(Jacobia{n)}^2$$

(39)

that is explained by the checkerboard ordering, thus the matrix, $\mathit{T}$, can be denoted as:

$$\left[\begin{array}{ll}0& {\mathit{T}}_1\\ {\mathit{T}}_2& 0\end{array}\right]$$

(40)

${\mathit{T}}_1$updates entries with odd entries, and ${\mathit{T}}_2$ updates odd entries with even entries. However, Successive Over-Relaxation significantly speeds up Gaussian-Siedel to solve Poisson-Boltzmann equation[80], and the optimal value to the relaxation parameter, $\omega =1/(1-\alpha )$:

$$\omega =2/(1+\sqrt{1-{\lambda }_N}$$

(41)

and the minimum spectral radius is:

$${\lambda }_N=({\lambda }_N^{gs}/{\left[1+\sqrt{1-{\lambda }_N^{gs}}\right]}^2)$$

(42)

${\lambda }_N^{gs}$denotes spectral radius of Gauss-Seidel. With the substitution, $\epsilon =1-{\omega }_N^{gs}$, the spectral radius:

$${\lambda }_N=(1-\epsilon )/{\left(1+\sqrt{\epsilon }\right)}^2 \approx 1+2\epsilon -2\sqrt{\epsilon }$$

(43)

and since $\epsilon$ is a small number, $\sqrt{\epsilon }$ is significantly larger. That denotes to significantly increase in rate of convergence. However, the spectral radius is sensitive to $\omega$, and the ${\lambda }_N$is chosen poorly, thus all improvements are lost [80], therefore Honig et. al. introduced an optimal Successive Over-Relaxation [80], which reduces the number of iterations by several orders of magnitude to a few iterations. However, ${\lambda }_N$​ from a few iterations, based on eigenvalue problems, finding the maximum eigenvalue $\mathit{T}$, is vital. This can be obtained by acquiring the minimum eigenvalue of $1-\mathit{T}$. That is related to Connected-Moments Expansion, which gives an approximation to the ground state wavefunction [80]. After approximation of the ground state wavefunction and Hamiltonian, and vector $\mathit{Q}$, is set to be zero:

$$\Phi (n)={\mathit{T}}^{\mathit{n}}\Phi (0)$$

(44)

Thus the prescription is to set the initial vector equal to the approximation highest eigenvalue eigenstate, set all charges equal to zero, and calculate a few moments [80].

Later, in order to better expression for Electrostatic Free Energy, they implemented a more accurate definition in DelPhI [81,82], which means that free energies became independent of the lattice used to solve the PB equation. The energy density can be expressed as:

$$∆g={\displaystyle \frac{1}{2}}{\rho }_{fixed}\Phi -{\sum }_i∆{\pi }_i-{\displaystyle \frac{1}{2}} {\rho }_{solv}\Phi$$

(45)

Additionally, in the finite-difference method, the system is discretized. The term ${\displaystyle \frac{1}{2}}{\rho }_{fixed}\Phi$, can also be expressed as:

$${\displaystyle \frac{1}{2}}{\rho }_{fixed}\Phi ={\displaystyle \frac{1}{2}}{\sum }_j{q}_j\rho \left({\overrightarrow{r}}_j\right)$$

(46)

where the electrostatic potential is calculated from all sources of charge except the one positioned at ${\overrightarrow{r}}_j$ [81,82]. This method eliminates accuracy concerns by representing the potential at a specific point as a superposition of real charges, surface polarization effects, and contributions from mobile ions, therefore, the potential at the position of charge $j$ can be written as:

$$\Phi \left({\overrightarrow{r}}_j\right)={\Phi }_{coul}\left({\overrightarrow{r}}_j\right)+{\Phi }_{reaction}\left({\overrightarrow{r}}_j\right)+{\Phi }_{ion}\left({\overrightarrow{r}}_j\right)$$

(47)

where Coulombic potential, generated by other fixed charges, can be expressed as:

$${\Phi }_{coul}\left({\overrightarrow{r}}_j\right)={\sum }_{i\ne j}{\displaystyle \frac{q_i}{4\pi {ϵ}_0{ϵ}_i\left({\overrightarrow{r}}_j-{\overrightarrow{r}}_j\right)}}$$

(48)

A corrected reaction field term, arising from the polarization of the boundary between different media is:

$${\Phi }_{reaction}\left({\overrightarrow{r}}_j\right)={\sum }_{i\ne j}{\displaystyle \frac{{\delta }_p}{4\pi {ϵ}_0{ϵ}_i\left({\overrightarrow{r}}_j-{\overrightarrow{r}}_j\right)}}$$

(49)

where $p$ denotes the positions of surface polarization charges, and ${\delta }_p$ represents the magnitude of the polarization charge at location $r_p$ on the surface. These charges are determined based on Gauss’s law and the electrostatic potential map generated by the finite difference method [81,82]. Last term, potential generated by mobile ions in solution:

$${\Phi }_{ion}\left({\overrightarrow{r}}_j\right)={\sum }_k{\displaystyle \frac{h^3{\delta }_{p_k^{ion}}}{4\pi {ϵ}_0{ϵ}_{solv}\left(r_j-r_k\right)}}$$

(50)

where $h$ denotes grid spacing, ${\delta }_{p_k^{ion}}$ denotes the net ion charge density at grid point $k$ in solution, which $k$ denotes grid points in the solution [81,82].

6. The APBS Model

The solution of Poisson-Boltzmann equation for biomolecules however remained as a challenge. Since the dielectric value of dielectric permittivity jumps drastically at the interface between solvent and solution; the function of the ionic strength of the solvent is discontinuous at the surrounding region of the biomolecule, and the formula itself features delta functions, and rapid nonlinearity that causes difficulties in analytical approaches. Different approaches have been proposed to approximate and calculate Poisson-Boltzmann Equation [84,85,86] in which some of them include uniform-mesh finite approaches [84,87], integral based methods [88,89], and nonadaptive finite element methods [86,90].

Here, we show a slightly different version of Poisson-Boltzmann equation for the 1:1 electrolyte:

$$-\nabla \left(ϵ\left(x\right)\nabla u\left(x\right)\right)+{\overline{K}}^2\left(x\right)\sinh \left(u\left(x\right)\right) ={\displaystyle \frac{4\pi e_c^2}{k_{B}T}}{\sum }_{i=1}^{N_m}{z}_i\delta (x-x_i)$$

(51)

where ${\overline{K}}^2$ is Debye–Hückel parameter, $N_m$ is the total number of point charges, $\delta (x-x_i)$ is the delta function that reflects point-charge behaivor of the charge at $x_i$ [91,92].

To address the challenges, Holst et. al. developed the Adaptive Poisson-Boltzmann Solver (APBS) model, where they adopted an adaptive multilevel finite element method to solve the nonlinear Poisson-Boltzmann equation, in which a posteriori error estimation guides the refinement [91,92]. Above written Poisson-Boltzmann equation is valid if the solution function differentiable twice, which is generally not satisfied in biomolecular systems. Therefore a weak formulation is required to the finite element method, due to obtaining reduced differentiability [91,92], that is:

$$〈F\left(u\right),v〉={\int }_{\mathsf{\Omega }}(ϵ\nabla u)\nabla v+{\overline{K}}^2\left(x\right)\sinh \left(u\right)v-fv dx$$

(52)

where, $v$ is the test function from Sobolev Space, $H_0^1(\Omega )$, and $f\left(x\right)$ represents the charge distribution term. In order to employ the Galerkin method, an approximation ${\overline{u}}_h$ to the boundary potential function$\overline{u}$ is constructed, and to satisfy Poisson-Boltzmann Equation [91,92]:

$$\mathrm{find} u_h\in {\overline{u}}_h+V_h\mathrm{such} \mathrm{that} F\left(u_h\right), v_h=0, \forall v_h\in V_h$$

(53)

where, $V_h$ is subspace of the Sobolev Space. Additionally, the Poisson-Boltzmann Equation follows as:

$$-\nabla \left(a\left(x\right)\nabla u\left(x\right)\right)+b\left(\mathrm{x},u\left(x\right)\right)=0 in \mathsf{\Omega },$$

(54)

$$u(x)=g(x) on \partial \mathsf{\Omega },$$

(55)

where, $a$ denotes $a:\mathsf{\Omega }⟼{\mathbb{R}}^{3\mathrm{x}3},\mathrm{b}:\mathsf{\Omega }\mathrm{x}\mathbb{R}⟼\mathbb{R},\mathrm{a}\mathrm{n}\mathrm{d} \mathrm{g}:\mathsf{\Omega } ⟼\partial \mathbb{R}$ [92].

In order to apply a Newton iteration, linearization is required [92]. To linearize the Poisson-Boltzmann equation operator,$〈F\left(u\right),v〉$, a bilinear linearization form $〈DF\left(u\right)w,v〉$ is produced:

$${\left.〈DF\left(u\right)w,v〉={\displaystyle \frac{d}{dt}}\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{m}〈DF\left(u+tw\right),v〉\right|}_{t=0}$$

(56)

and it follows,

$$={\int }_{\mathsf{\Omega }}a\nabla w\nabla v+{\displaystyle \frac{\partial b(x,u)}{\partial u}}vw dx$$

(57)

where $u$, $v$, and $w$are three arguments of the scalar-valued function. $〈F\left(u\right),v〉$, and the associated bilinear linearization form, $〈DF\left(u\right)w,v〉$ together with a continuous piecewise polynomial subspace of the solution space, ${\overline{u}+H}_0^1(\Omega )$, are required to employ the finite element method for numerical solution of the original elliptic equation [92].

Moreover, the inexact Newton method provides a numerical solution to the nonlinear Poisson-Boltzmann equation by iteratively solving its linearized versions to calculate adjustments until a desired error tolerance is met, where computationally intensive part of this method is complexity of the linear algebraic equations generated in each step. Notably multilevel methods are near-optimal techniques to solving such equations resulting from discretization of a large class of general linear elliptic problems [92]. These methods operate by using an interpolation operator to map the discretized linearized equation onto a coarser mesh, which reduces the number of variables and simplifies the problem. The approximate solution on the finest mesh can then be obtained by using the coarse solutions to accelerate the convergence of iterative methods on the finer levels [91,92].

Other aspect of the APBS is the posteriori error estimation and adaptive refinement, where the main idea can be understood by considering an approximation, $\tilde{u}$, to the solution of a linear equation defined by a nonsingular operator A [91]:

$$Au=f$$

(58)

The quality of the approximation by forming the residual can be tested and since error in the approximation $e = u- \tilde{u}$ satisfies the error equation:

$$Ae=A\left(u-\tilde{u}\right)=Au-A\tilde{u}=f-A \tilde{u}=r$$

(59)

Inverting A and taking norms of both sides, thus gives:

$$\left|\left|e\right|\right|=\left|\left|{\mathrm{A}}^{-1}r\right|\right|\le \left|\right|A^{-1}\left|\right|\left|\right|\mathrm{r}\left|\right|$$

(60)

This inequality shows that the error in approximation is limited by the product of the inverse operator and residual norm. Thus, estimation of error of the solution can be obtained [91,92].

7. The ABSINTH Model

It has been noted that, due to EEF1’s incapabilities, other implicit solvent approaches have advantages such as in reproducing the disordered nature of an amyloid fragment [93]. Intrinsically disordered proteins are widely present in the archaeaic, prokaryotic and eukaryotic proteome [94,95,96]. Due to the fact that implicit solvents are prone to fail to model these disordered structures, Pappu et. al. introduced the self-Assembly of Biomolecules Studied by an Implicit, Novel, and Tunable Hamiltonian (ABSINTH) model [28]. Although it is related to the work of Lazardis et. al., it was developed to address accurate modeling of intrinsically disordered proteins. The Hamiltonian of this model can be formulated as:

$$E_{total}=W_{Solv}+U_{LJ}+W_{ele}+U_{corr}$$

(61)

$W_{Solv}$is the solvation term that corresponds to the direct mean field interaction, $U_{LJ}$ represents the dispersive and short-range steric interactions contributions, $U_{corr}$ defines the torsional correction terms applied to dihedral angles influenced by electronic effects which couldn’t be captured by $U_{LJ}$, and $W_{el}$ represents electrostatic interactions modulated by the mean-field dielectric [28]. Since the degree of the solvent accessibility has influence on the direct mean field interactions, it can be formulated as:

$$W_{solv}={\sum }_{i=i}^{N_{SG}}{\zeta }_i∆G_{solv}^i={\sum }_{i=i}^{N_{SG}}\left[{\sum }_{k=1}^{n_i}{\lambda }_k^i{v}_k^i\right]∆G_{solv}^i$$

(62)

Accordingly, it denotes the solvation free energy of solvation group i, while the average weighted energy over its possible states, $v_k^i$, is calculated using the weight factors weight factors ${\lambda }_k^i$ [28].

8. Quasi-Chemical Theory

Quasi-Chemical theory can be summarized as a molecular level framework for estimating the excess chemical potential of solutes, where the region around the solute of interest is partitioned into inner- and outer-shell domains [97,98,99]. For instance, it implies that the chemical potential of a lithium ion ${(\mu }_{Li}^{+})$can be expressed in terms of interaction of ideal and non-ideal contributions [98]:

$${\beta }^{*}{\mu }_{L{i}^{+}}= {\log }_e({\displaystyle \frac{{\rho }_{L{i}^{+}}}{q_{L{i}^{+}}/V}})-{\mathrm{l}\mathrm{o}\mathrm{g}}_e{\rho }_0- {\mathrm{l}\mathrm{o}\mathrm{g}}_e({\sum }_{n\ge 0}{Ǩ}_n\rho H_2{O}^n)$$

(63)

We should note here that ${\beta }^{*}$ differs from $\beta$ as it is equal to ${\displaystyle \frac{1}{RT}}$ unlike β.

The ion hydration reactions can be modeled through inner shell interactions without considering medium effects [98,100]:

$$L{i}^{+}+nH2O \rightleftharpoons L{i}^{+}{\left(H_{2}O\right)}_n$$

In order to understand the Quasi-Chemical Theory, we need to describe the potential distribution theorem [98,101]:

$${\rho }_{\sigma }={⟨e^{-{\beta }^{*}∆U}⟩}_0{\mathcal{z}}_{\sigma }(q_{\sigma }/V)$$

(64)

${\rho }_{\sigma }$ is the density of the considered molecule, $q_{\sigma }$ is the single molecule partition function for that molecule, ${\mathcal{z}}_{\sigma }$ is the absolute activity of that molecule and equals to;

$${\mathcal{z}}_{\sigma }=e^{{\beta }^{*}{\mu }_{\sigma }}$$

Asthagiri et. al studied the influence of chemical effects which are based on inner shell reactions [97]:

$$X^{\pm }+n{H}_{2}O \rightleftharpoons X{{\left(H_{2}O\right)}^{\pm }}_n$$

$X^{\pm }$ denotes various ions such as $H^{+}, L{i}^{+}, N{a}^{+}, O{H}^{-}$ [97].

Furthermore, quasi-chemical approximation to the excess chemical potential of an ion in water can be described as:

$$\beta *{\overline{\mu }}_{X\pm \left(w\right) }\approx -\ln ({Ǩ}_n{\rho }_{H_2{O}^n})$$

(65)

and ${Ǩ}_n$ equals to:

$${Ǩ}_n= {Ǩ}_n^{(0)}*e^{-\beta \left({\overline{\mu }}_{X{\left(H_{2}O\right)}_n}^{\pm }\left(w\right)-{\overline{\mu }}_{H_{2}O}(w)\right)}$$

${Ǩ}_n^{(0)}$refers to the equilibrium constant for the reaction in an ideal gas state, and ${\overline{\mu }}_{X{\left(H_{2}O\right)}_n}^{\pm }\left(w\right)$ is the excess chemical potential of the ion–water cluster, and ${\overline{\mu }}_{H_{2}O}(w)$ is the excess chemical potential of a water molecule in bulk water.

9. Transfer Free Energy Approach

(Mis)folding of intrinsically disordered proteins remains a challenge [102]. Shea et. al. introduced more accurate implicit solvation model for such systems [29].

The model can be defined as:

$$G_{solv}=E_{vac}+G_{ele}+G_{np}$$

(66)

$E_{vac}$ represents the peptide energy in vacuum which results from both internal bonded contributions such as angles and dihedrals and non-bonded van der Waals interactions. They also introduced a new approach where the surface tension coefficient has a dependency on temperature. In fact, the probability of not buried side chain of residue $i$, $P_{nb}^{sc,i}$, can be calculated as:

$$P_{nb}^{sc,i}={\displaystyle \frac{1}{N_i}}{\sum }_{k=1}^{N_i}\left[{\displaystyle \frac{SA{S}_k^{sc,i} }{SA{S}_{max}^{sc,i} }}\right]$$

(67)

$SA{S}_k^{sc,i}$ defines side chain solvent accessibility of the $k^{\text{'}}th$ residue type $i$ and $SA{S}_{max}^{sc,i}$ refers to the maximum solvent accessibility of side chain type i [30] Based on this logic. The probability for the backbone can be calculated as:

$$P_{nb}^{bb}={\displaystyle \frac{1}{N_i}}{\sum }_{k=1}^{N_i}\left[{\displaystyle \frac{SA{S}_k^{bb} }{SA{S}_{max}^{bb} }}\right]$$

(68)

Thus, the energy can be calculated as:

$$E_{nb}^{sc,i/bb }=-RTln{P}_{nb}^{sc,i/bb }$$

(69)

The buried side chain can be calculated using:

$$E_b^{sc,i/bb }=-RTln(1-P_{nb}^{sc,i/bb })$$

(70)

These can help in calculating the energetic cost for exposing side chain type $i$ or the peptide backbone between two different temperatures via:

$${\left(E_{nb}^{sc,i/bb}-E_{bb}^{sc,i/bb}\right)}_T- {\left(E_{nb}^{sc,i/bb}-E_{bb}^{sc,i/bb}\right)}_{T_0}=[{\gamma }_2^{sc,i/bb}(T)-{\gamma }_0]SA{S}_k^{sc,i/bb}$$

(71)

Nonpolar contribution can be expressed as:

$$G_{np}=\left\{\begin{array}{l}{\gamma }_{0}SAS, T=T_0\\ {\gamma }_{0}SAS+{\sum }_{k=1}^n({\gamma }_2^{sc,k}(T)-{\gamma }_0)SA{S}^{sc,k} +{\sum }_{k=1}^n\left({\gamma }_2^{bb}\left(T\right)-{\gamma }_0\right)SA{S}^{bb,k} , T\ne T_0 \end{array}\right.$$

(72)

Fractional solvent accessibility can be described as;

$$G_{np}=\left\{\begin{array}{l}{\gamma }_{0}SAS, T=T_0\\ {\gamma }_{0}SAS+{\sum }_{k=1}^n∆g_{tr,2}^{sc,k}\left(T\right){\alpha }^{sc,k}+∆g_{tr,2}^{bb}\left(T\right){\sum }_{k=1}^n{\alpha }^{bb,k} , T\ne T_0 \end{array}\right.$$

(73)

where, ${\alpha }^{sc,k/bb,k}$ is a fractional solvent accessibility parameter that can be found in Ref.30 and $∆g_{tr,2}^{sc,k/bb,k}$ can be written as:

$$∆g_{tr,2}^{sc,k/bb,k}=({\gamma }_2^{sc,k/bb}(T)-{\gamma }_0)SA{S}_{GLY-k-GLY}^{sc,k/bb,k}$$

The general formula can be written as:

$$G_{solv}=E_{vac}+G_{ele}+G_{np}^0+G_{tr}(T)$$

(74)

The second approach can be described as:

$$G_{tr}(T)=\left\{\begin{array}{l}0, T=T_0\\ {\gamma }_{0}SAS+{\sum }_{k=1}^n∆g_{tr,2}^{sc,k}\left(T\right){\alpha }^{sc,k}+ ∆g_{tr,2}^{bb}\left(T\right){\sum }_{k=1}^n{\alpha }^{bb,k} , T\ne T_0 \end{array}\right.$$

(75)

$${G_{np}^0=\gamma }_{0}SAS$$

(76)

Shea et. al. [29] also developed a recent approach based on the one described above, due to the need for accurate folding of both natively ordered proteins and intrinsically disordered proteins calculations [102]. Their work demonstrated that the approach, which we describe above, was prone to produce proper ensembles for disordered peptides [29]. Their third approach is defined as [29]:

$$G_{tr}(T)={\sum }_{k=1}^n∆g_{tr,2}^{sc,k}\left(T\right){\alpha }^{sc,k}+ ∆g_{tr,2}^{bb}\left(T\right){\sum }_{k=1}^n{\alpha }^{bb,k} , \forall T$$

(77)

$$G_{np}^0=0$$

(78)

This approach also failed in calculating proper results since it overestimates the helical propensity of disordered peptides. Thus, they demonstrated that a combined approach shows more promising results:

$${G_{np}^0=\gamma }_{0}SAS$$

(79)

10. The GBNSR6 Model

The GBNSR6 is a grid-based molecular surface implementation that uses the R6 variant of the generalized Born implicit model. The R-6 (R6) integral allows the calculation of Born radii. GBNSR6 in AmberTools is integrated into MMPBSA.py and can be used to calculate the solvation free energy [103].

Izadi, Aguilar, and Onufriev used the GBNSR6 model to predict the protein-ligand binding energies. Specifically, the binding free energies of 15 small protein-ligand complexes obtained using GBNSR6 were examined by utilizing TIP3P and TIP4Ew OPC explicit solvent free energies as reference. The root mean square deviation of GBNSR6 using TIP3P was found to be close to the error margin of explicit solvents. $\mathsf{\Delta }\mathsf{\Delta }{G}_{pol}$calculated with GBNSR6 was closer to the calculations utilizing the OPC model. It was shown that almost all deviations are improved by a simple uniform scaling of the set of Bondi radii. Large discrepancies were found between the binding and solvation energies using explicit models. Thus, the need for developing more accurate implicit solvent models was emphasized. We should note that it was stated that it would be more accurate to measure the model accuracy by comparing with experimental data [104].

In a study by Wang et al., the performance of various MM/PB(GB)SA approaches for predicting protein–protein interactions and protein–ligand binding structures was systematically evaluated using a dataset of 900 unique docking poses. The MM/PBSA method demonstrated superior pose ranking accuracy compared to MM/GBSA and Glide SP scoring, with the PB3 model achieving a success rate of 74%. Among the MM/GBSA approaches, the GB6 model was particularly notable, highlighting the effectiveness of the GBNSR6 implicit solvent model. Additionally, this work led to the development of the Fast Amber Rescoring for PPI Inhibitors (farPPI) web server, which enables efficient scoring of docking poses using MM/PB(GB)SA methods [105]. In the study by Forouzesh et al., the impact of different atomic radii parameterizations on binding free energy calculations was assessed using the Ras-Raf protein complex, and the resulting absolute binding free energies showed good agreement with experimental values. Subsequently, the performance of the GBNSR6 model was evaluated for predicting the binding free energy between the SARS-CoV-2 spike receptor binding domain (S RBD) and the human ACE2 receptor. When using the Bondi radii set, the calculated binding free energy was overestimated, while the OPT1 radii set led to an underestimation. For the SARS-CoV-2 S RBD–ACE2 complex, the computed ΔG_bind values deviated by approximately 4 kcal/mol from the experimental reference, indicating limited quantitative accuracy. Notably, the two ΔG_bind estimates for this complex suggest a near cancellation between the large enthalpic (ΔH) and entropic (−TΔS) contributions to binding free energy, supporting the qualitative reliability of the MMGB/SA approach in such cases [106].

Moreover, Cain, Risheh and Forouzesh presented a Physics-Guided Neural Network (PGNN) model for protein-ligand binding (despite the possibility of overfitting of data-driven models). This model takes its physics basis from GBNSR6 and its data model from Graph Convolutional Network. This physics-based AI model predicted the binding free energy more accurately than GraphConv and AtomicConv data-driven models. It showed consistent results in both training and test sets, demonstrating its transferability between different data sets. In order to measure the robustness of the PGNN model, a noise function was added and evaluated. Since adding or not adding noise did not significantly change the results, the robustness of the model was confirmed. When looking at the physical properties, it was seen that the van der Waals energy was not correlated with other properties. It was suggested that this property should be preserved and the other properties should be eliminated or grouped [107].

Furthemore, Tolokoh et al. studied the implicit water multipole GB (IWM-GB) model that was developed based on the GBNSR6. In addition to incorporating dipole water polarization as in GBNSR6, this model also accounts for the multipolar field effects of water molecules within the first hydration shell of the solvent. The model showed results compatible with experimental data. Root mean squared error decreased by 12% compared to the usage of the TIP3P model. The IWM-GB model was examined in two sub-versions as IWM-GB WC and IWM-GB NC. While IWM-GB WC includes the polar-nonpolar interaction, IWM-GB NC does not. The charge hydration asymmetry (CHA) effect of the WC model was noted to be more suitable for achieving experimental results. Thus, the importance of the polar-nonpolar interactions in the model was demonstrated. The IWM-GB model does not take into account the hydrogen bridges formed by water molecules between two atoms in polypeptides [108].

11. Implicit Quantum Models:

As mentioned above, the Born model was used to calculate the solvation energy of a spherical ion in a dielectric medium [109]. The model is based on classical electrostatics and remains valid at the point where the solvent is considered as a continuous dielectric medium [110]. The Born model offers a straightforward and intuitive formula for solvation energy, which is proportional to the square of the charge and inversely related to the ion's radius [111]. However, its use is limited to spherical ions. It also does not take into account the shape of molecules or non-electrostatic factors such as cavitation. The Poisson Boltzmann equation, considered an important model for electrostatic solvation, describes the electrostatic potential in a medium in which a charged solute is surrounded by a dielectric solution [109]. GB model, developed in the 1980s, builds upon the Born model to accommodate molecules of various shapes [36].

This is accomplished by introducing effective Born radii, which estimate the solvation energy of a molecule by summing pairwise interactions between atomic charges [9]. The GB model is known to be computationally efficient and appropriate for large systems such that it is quite a popular parameter model for molecular dynamics (MD) simulations [112]. The validity of this model is in doubt for some systems as it depends on the Born radii used and the pairwise interaction approximation [113]. The 1980s saw also an important step forward in implicit solvent model development. In polarizable continuum model (PCM), the solvent is treated as a polarizable continuum and the electrostatic solute-solvent interaction is estimated by solving the Poisson equation [114]. The PCM solvation energy is written as:

$$\mathsf{\Delta }{G}_{solv}={\displaystyle \frac{1}{2}}\int \rho \left(r\right)\varphi reaction\left(r\right)dr$$

(80)

Here, ρ(r) is the charge density of solute and $\varphi reaction\left(r\right)$ is the reaction potential due to solvent. The Poisson equation is solved to obtain the reaction potential:

$$\nabla \left(\epsilon \left(r\right)\nabla {\varphi }_{reaction}\left(r\right)\right)=-4\pi \rho \left(r\right)$$

(81)

where $\epsilon \left(r\right)$ is the solvent dielectric constant. PCM gives a more accurate description of the solute-solvent interaction than the Born and GB models, since the detailed charge distribution of the solute is taken into account [12]. However, the numerical resolution of the Poisson equation can become computationally costly, especially for large systems [115].

The Conductor-like Screening Model (COSMO) is a simplification of the PCM approach, which assumes that the solvent behaves as a conductor and has an infinite dielectric constant. This is based on the assumption that increases the computational efficiency of COSMO at the same time preserving a reasonable accuracy. Solvation energy in COSMO is described by [116]:

$$\mathsf{\Delta }{G}_{solv}={\displaystyle \frac{1}{2}}\sum {\sigma }_i{\varphi }_i$$

(82)

Here, ${\sigma }_i$ is the surface charge density at point i on the solute-solvent boundary, and ${\varphi }_i$ is the electrostatic potential at point i; the surface charges are given by the solution of the matrix equation [117].

$$A\sigma =B\varphi$$

(83)

COSMO is commonly used method in DFT calculations due to its computational ease [118]. It has an infinite dielectric constant however, and this is a perfect approximation for solvents in many cases but expands sources of error when using solvents with a lower dielectric constant (nonpolar ones) [119]. In contrast, for non-polar systems, C L-based COSMO overestimates the solvation energy, which leads to less accurate predictions [120]. The underlying reason for this limitation is the inability of COSMO to describe the weak electrostatic and specific solute-solvent interactions (hydrogen bonding, van der Waals) that exist in non-polar media. To circumvent these limitations, more comprehensive models such as COSMO-RS (Real Solvent) have been developed [13]. Extensions by statistical thermodynamics that take into account different solute/solvent interactions such as hydrogen bonding and van der Waals forces are built upon the original COSMO model, which lead to the COSMO methods (Continuum solvation COSMO) and the COSMO-RS [121]. Therefore, compared to QSPR, COSMO-RS has the potential to more accurately summarize predictions for a much broader range of solvents, such as ionic and non-polar liquids. COSMO-RS is a versatile method and can also estimate thermodynamic properties such as activity coefficients, vapor-liquid equilibria, and partition coefficients [122].

Due to its simplicity and computational efficiency, COSMO is still commonly used for solvation studies in a high-dielectric and polar solvents, despite some discrepancies. Performance and scalability are important for high throughputs application and DFT calculations where this is very useful. For higher accuracy, the use of more sophisticated models, like COSMO-RS or SMD (Solvation Model based on Density), is recommended in systems with non-homogeneous, non-polar solvents or complex solute-solvent interactions [123]. In fact, COSMO-RS uses statistical thermodynamics to account for some solute-solvent interactions, such as hydrogen bonding, van der Waals forces, and solvent-solvent interactions, whereas COSMO assumes the solvent acts as a conductor with an infinite dielectric constant [124]. For predicting thermodynamic parameters like activity coefficients, vapor-liquid equilibria, and solubility in a variety of solvents, including polar, non-polar, and ionic liquids, COSMO-RS is therefore more accurate [125]. In order to account for non-electrostatic interactions, COSMO-RS incorporates a residual term into the electrostatic framework of COSMO.⁸¹ There are two primary components to the overall solvation energy in COSMO-RS [126]; when considering the electrostatic component side, this component is computed using the same methodology as in COSMO, which makes use of the electrostatic potentials ${\sigma }_i$ and surface charge densities ${\varphi }_i$ [125]:

$$∆G_{COSMO}={\displaystyle \frac{1}{2}}\sum {\sigma }_i{\varphi }_i$$

(84)

The solvent's polarization effects are captured by this phrase. This component takes into consideration non-electrostatic interactions on the residual element side, including solvent-solvent interactions, van der Waals forces, and hydrogen bonding. The expression for the remaining energy is [125].

$$∆G_{res}={\sum }_{ i,j}{\mathsf{\Gamma }}_{i,j}{\sigma }_i{\sigma }_j$$

(85)

In this case, ${\sigma }_i$ and ${\sigma }_j$ are the surface charge densities at surface segments i and j, and ${\mathsf{\Gamma }}_{ij}$ is the interaction energy between these segments [127]. COSMO-RS, model has a number of noteworthy benefits that make it an effective tool for solvation modeling and thermodynamic property prediction [124]. Its extensive application to a variety of solvents, including polar solvents like alcohols and water, non-polar solvents like hexane and toluene, and even ionic liquids, is one of its main advantages [124].

In contrast to the original COSMO model, which is mostly accurate for high-dielectric solvents and assumes an infinite dielectric constant, COSMO-RS takes into account particular solute-solvent interactions like solvent-solvent interactions, hydrogen bonding, and van der Waals forces [128]. This allows it to produce highly accurate predictions for systems like complex mixes or non-polar solvents where these interactions are crucial [129]. The capacity of COSMO-RS to accurately forecast thermodynamic characteristics is another significant benefit [130]. It is extremely useful in drug discovery and (bio)chemical process design since it is frequently used to compute activity coefficients, vapor-liquid equilibria, solubility, and phase behavior [124]. COSMO-RS also achieves a balance between computing efficiency and accuracy [131]. The addition of residual factors for non-electrostatic interactions makes it more computationally demanding than COSMO, but it is still effective enough to manage large systems. For comparison purposes, traditional COSMO and COSMO-RS can be compared as follows [132].

The SMx models, which were developed from the 1990s to the 2000s, represent a family of semi-empirical solvation models that incorporate both electrostatic and non-electrostatic contributions [33]. These models are designed to be universal, allowing them to be applied across a variety of solvents and solutes [33]. The solvation energy in the SMx models is expressed as [33]:

$$∆G_{solv}=∆G_{elec}+∆G_{non-elec}$$

(86)

or,

$$G_{solv}=G_{elec}+G_{cav}+G_{disp}+G_{rep}$$

(87)

The electrostatic term $∆G_{elec}$is calculated using a dielectric continuum model like PCM or GB, while the non-electrostatic term $∆G_{non-elec}$ encompasses cavitation, dispersion, and repulsion effects. The non-electrostatic term is formulated as [33]:

$$∆G_{non-elec}=\gamma A+\alpha B+\beta$$

(88)

γ, α, and β are parameters, while A is the solvent accessible surface area. SMx models are frequently employed in quantum modeling and are renowned for their exceptional precision. However, the quality of the parameterization determines how accurate they are [33]. M5.4, SM6, SM8, and SMD are among the models in the SMx family; they differ in their mathematical formulations and parameterizations [133]. To study certain solvation modeling issues, these variations are required [134]. The balance between efficiency and precision is an important factor [135]. Models like SM5.4 and SM6 are perfect for high-throughput screening of large molecular datasets because they prioritize computational efficiency by using simplified equations for non-electrostatic variables [136]. However, by incorporating more thorough parameterization and empirical terms—both of which are essential for precisely computing solvation free energy in complicated systems—models such as SM8 and SMD place an emphasis on accuracy [137]. Their compatibility for various solvents and solutes is another crucial factor.

Certain models are optimized for specific solutes, like ions or neutral molecules, or for specific solvent types, like water or organic solvents [33]. For instance, SM5.4 is especially made for tasks like analyzing neutral molecules in aqueous solution, whereas SMD is a general-purpose model that can handle a wide range of solvents and solutes, including ionic species [33]. Additionally, choosing a model is influenced by compatibility with quantum mechanical techniques [33]. As a matter of fact, the SMx models may be used to DFT and other quantum mechanical techniques, including semi-empirical techniques such as AM1 and PM3 [136]. Model selection is impacted by the theoretical level used in the computations [15]. For instance, SMD and DFT are commonly used for high-accuracy calculations, while SM5.4 and semi-empirical methods are often used for faster calculations [137]. Because of these variations in design and use, the SMx family is flexible and able to handle a variety of computational chemistry problems [15].

Table 1. A summary of SMx models and their characteristics.

Model	Key Feature	Mathematical Approach	Advantages	Applications
SM 5.4	Designed for fast and efficient calculations.	Uses the GB model for electrostatic terms. Non-electrostatic terms (cavitation, dispersion, repulsion) are modeled using surface area and volume.	Computationally efficient. Suitable for high-throughput screening.	Studying neutral molecules in aqueous solutions.
SM 6	An improved version of SM5.4. Includes more accurate parameterization	Uses PCM or GB for electrostatic terms. Non-electrostatic terms are modeled with more detailed parameters.	Higher accuracy than SM5.4. Still computationally efficient.	Studying ionic solutions and solubility of organic molecules.
SM 7	Focuses on improving non-electrostatic terms.	Uses PCM for electrostatic terms. Non-electrostatic terms are modeled with advanced empirical parameters.	Better accuracy for non-electrostatic effects.	Solvation of polar and non-polar molecules in various solvents.
SM 8	A universal model optimized for both ionic and neutral molecules.	Uses PCM for electrostatic terms. Non-electrostatic terms are modeled using surface area, volume and advanced empirical parameters.	High accuracy for a wide range of solvents and solutes.	Drug design, ionic solutions, and solubility of organic molecules.
SMD	The most advanced SMx model. Designed as a universal solvation model.	Uses PCM for electrostatic terms. Non-electrostatic terms are modeled using surface area, volume and empirical parameters (γγ, αα, ββ).	High accuracy across a wide range of solvents and solutes. Compatible with DFT.	Drug design, materials science, environmental chemistry (e.g., solubility of pollutants).
SM 12	An extension of SMD with improved parameterization	Uses PCM for electrostatic terms. Non-electrostatic terms are modeled with more refined empirical parameters.	Enhanced accuracy for specific solvent-solute systems.	High-precision calculations for solvation free energies in complex systems.
SM x-NP	Designed for non-polar solvents and solutes.	Uses GB or PCM for electrostatic terms. Non-electrostatic terms are modeled with parameters optimized for non-polar interactions.	Accurate for non-polar systems.	Studying organic semiconductors and polymers in non-polar solvents.
SM x-IL	Tailored for ionic liquids.	Uses PCM for electrostatic terms. Non-electrostatic terms are modeled with parameters optimized for ionic liquids.	High accuracy for ionic liquid systems.	Studying solvation and reactivity in ionic liquids

Table 1. A summary of SMx models and their characteristics.

Model	Key Feature	Mathematical Approach	Advantages	Applications
SM 5.4	Designed for fast and efficient calculations.	Uses the GB model for electrostatic terms. Non-electrostatic terms (cavitation, dispersion, repulsion) are modeled using surface area and volume.	Computationally efficient. Suitable for high-throughput screening.	Studying neutral molecules in aqueous solutions.
SM 6	An improved version of SM5.4. Includes more accurate parameterization	Uses PCM or GB for electrostatic terms. Non-electrostatic terms are modeled with more detailed parameters.	Higher accuracy than SM5.4. Still computationally efficient.	Studying ionic solutions and solubility of organic molecules.
SM 7	Focuses on improving non-electrostatic terms.	Uses PCM for electrostatic terms. Non-electrostatic terms are modeled with advanced empirical parameters.	Better accuracy for non-electrostatic effects.	Solvation of polar and non-polar molecules in various solvents.
SM 8	A universal model optimized for both ionic and neutral molecules.	Uses PCM for electrostatic terms. Non-electrostatic terms are modeled using surface area, volume and advanced empirical parameters.	High accuracy for a wide range of solvents and solutes.	Drug design, ionic solutions, and solubility of organic molecules.
SMD	The most advanced SMx model. Designed as a universal solvation model.	Uses PCM for electrostatic terms. Non-electrostatic terms are modeled using surface area, volume and empirical parameters (γγ, αα, ββ).	High accuracy across a wide range of solvents and solutes. Compatible with DFT.	Drug design, materials science, environmental chemistry (e.g., solubility of pollutants).
SM 12	An extension of SMD with improved parameterization	Uses PCM for electrostatic terms. Non-electrostatic terms are modeled with more refined empirical parameters.	Enhanced accuracy for specific solvent-solute systems.	High-precision calculations for solvation free energies in complex systems.
SM x-NP	Designed for non-polar solvents and solutes.	Uses GB or PCM for electrostatic terms. Non-electrostatic terms are modeled with parameters optimized for non-polar interactions.	Accurate for non-polar systems.	Studying organic semiconductors and polymers in non-polar solvents.
SM x-IL	Tailored for ionic liquids.	Uses PCM for electrostatic terms. Non-electrostatic terms are modeled with parameters optimized for ionic liquids.	High accuracy for ionic liquid systems.	Studying solvation and reactivity in ionic liquids

The Solvation Model based on Density (SMD) is a model that merges the advantages of PCM for electrostatics with empirical non-electrostatic terms. [123] SMD is intended for use with a broad spectrum of solvents and solutes, including both neutral molecules and ions. The solvation energy in SMD is represented as;

$$∆G_{solv}=∆G_{elec}+\gamma A+\alpha B+\beta$$

(89)

The electrostatic term $∆G_{elec}$ is derived from a dielectric continuum model, while the non-electrostatic terms γA+αB+β account for cavitation, dispersion, and repulsion effects. SMD is widely used to forecast solvation free energies and is compatible with a number of quantum mechanical techniques, such as DFT [139,140]. However, the parameterization for various solvents affects its accuracy, just like it does for other SMx models [14]. The quality of the dataset used for parameterization determines the accuracy of the empirical parameters ($\gamma ,\alpha and \beta$), which are obtained from experimental data [33]. Because the settings have been well adjusted using a sizable database of experimental solvation free energies, SMD is extremely accurate for solvents like water and organic solvents [33]. However, the parameterization might be less accurate for ionic liquids or less widely used solvents. SMD uses a parameterization technique that attempts to cover a wide range of solvents and solutes in order to overcome this restriction. With this method, a variety of experimental data, such as solvation free energies, activity coefficients, and partition coefficients, are fitted to empirical parameters [33].

We should note that the solvent-solute system can still affect SMD accuracy [141]. For example, certain solute-solvent interactions, such as hydrogen bonding or van der Waals forces, may be difficult for the model to describe in ionic liquids or strongly non-polar solvents [142]. Apart from the difficulties in parameterization, the accuracy of SMD is also affected by the chosen quantum mechanical method [140]. Although DFT and semi-empirical approaches are compatible with SMD, the quality of the electrostatic and non-electrostatic terms may be impacted by the level [143]. High-level DFT simulations with precise electron density distributions, for example, typically produce better results than semi-empirical approaches, which might incorporate more approximations [143]. Despite these difficulties, SMD's adaptability and applications make it one of the most utilized models. It is widely used in drug design, where solvation free energy estimates are essential for determining binding affinities and solubility.

12. Some Applications in Biology

Intrinsically disordered proteins are proteins that do not have a stable three-dimensional structure [144]. Techniques such as X-ray crystallography and cryo-electron microscopy are insufficient in explaining their disordered nature. Their properties are investigated with applications combined with implicit solvent models. Different dissolution times are observed by representing continuous products instead of individual atoms. When working with existing implicit solvent models; the formation of excessively collapsed disordered states in explicit solvent models has been observed because of comparisons made with FRET and SAXS experiments. The problem of overestimating the α-helix and β-sheet structures has been largely solved by the change of the backbone torsional change based on the new generation NMR systems [145]. For instance, the EEF1 model by Karplus and Lazaridis uses the solvation free energies in empirical runs adjusted according to the degree of buriedness of the protein functional groups. It calculates the electrostatic interactions in solution with the distance dependent dielectric constant by neutralizing the ionic side chains. The EEF1 model separates the free energy of a protein in solution into two main components: intramolecular energy and solvation free energy. The interaction of an atom with solvent is measured by how much the surrounding atoms restrict solvent access [79]. The ABSINTH model by Vitalis and Pappu models these proteins better; EEF1 and ABSINTH differ in the choice of solvation groups and the way they measure solvent accessibility. This model has a two-component approach to describe the transfer of solute to the continuum, modeling the DMFI (direct mean-field interaction) and the screening of polar interactions. The polar and nonpolar parts of the transfer process are handled using reference solvation free energies for the solvation groups [28].

Accurate calculation of protein-ligand binding energy is an essential step in drug design. There are many implicit solvent models used in this field. Some of them are PCM, S-GB (Surface Generalized Born), COSMO , GBNSR6 (Generalized Born R6 version) and PB model [146]. In a study conducted by Katkova et al., hydration, solvation and desolvation energies calculated with these models were compared with the results calculated using explicit solvent models and experimental data. They focused on polar components rather than non-polar components. As a result, although Poisson-Boltzmann and GBNSR6 achieved the most accurate results in terms of accuracy in calculating desolvation energies, GBNSR6 model stood out in terms of speed. However, it was seen that better parameter selections are still needed in desolvation energy calculations. This is because the choice of atomic radii and charges has a strong influence on the calculation of desolvation energy. None of these models could provide chemical accuracy (error level less than 1 kcal/mol) [146].

Guo, Zuojun, et al., used Coloumb area approximation (CFA) and Level-Set Variational Implicit-Solvent Model (VISM) for determining protein binding surfaces. Most small molecules tend to bind to hydrophobic pockets with complex topology on the surface of proteins. In this method, the dissolution process of molecules was described by minimizing the solvation free energy of solute-solvent interfaces. Surface tension, electrostatic and van der Waals interactions were calculated using an implicit solvent model. In this way, the depth, volume and hydrophobic properties of binding pockets could be rapidly determined [147]. Furthermore, Feig and co-workers performed DNA and DNA-protein simulations using the GBMV model. The effect of salt presence was also evaluated. Obtained results were analyzed by comparing them with explicit solvent simulations and experimental data. Simulations performed with the default radius set in CHARMM were calculated using root mean square deviations. The DNA structure remained stable throughout the simulation. The simulated structure was quite similar to the experimental structure. When the standard helical parameters were examined, it was seen that they are in good agreement with the average helical parameters obtained from explicit solvent simulations, but the standard deviations are two times larger. Since implicit solvent models capture the effects of water molecules in less detail, a higher degree of oscillation was observed [148].

In a study by Prabhu et al., a hybrid implicit water solvation model was developed for DNA and RNA molecular simulations. Specifically, Finite Difference Poisson–Boltzmann (FDPB) and GB models were generally applied in cases where molecules with small net charges and ionic strength effects are not dominant. This hybrid model was developed due to the lack of detailed examination of the cases where molecules carry high charges. In the FDBP model, the ionic strength was set to zero, thus reducing it to the finite-difference Poisson (PDB) model. Ions were explicitly treated as solvents. Since ions constitute a smaller part of the total number of the system, this explicit approach did not increase the computational power significantly. It was observed that long-term stability was provided in simulations performed on B-DNA dodecamers and decamers and tetraloop RNA molecules. Similar results were obtained for explicit water simulations. In the analyses performed with the CURVES DNA structure analysis program, it was observed that there was no statistically significant difference in 24 helicoidal and base parameters. The implicit solvent model accurately simulated the transition from A to B, and similar final stable B structures were achieved with previous studies. Results were obtained faster compared to the usage of explicit solvent models [149].

13. Future Perspectives

Implicit solvent models have been useful in biomolecular simulations and quantum mechanical calculations, yet ongoing advances in methodologies and our understanding of solvation phenomena promise to drive the field forward. A key future direction involves the integration of artificial intelligence to refine parameterization, improve transferability across diverse environments, and accelerate the development of next-generation hybrid models. Physics-guided neural networks and data-driven approaches can further enhance the predictive accuracy and computational efficiency of implicit solvation methods. Quantum computing and multi-scale modeling frameworks offer exciting opportunities for simulating large, complex biomolecular systems with accuracy. Combining quantum models with high-level electronic structure methods will enable the study of challenging processes such as enzymatic catalysis and the behavior of intrinsically disordered proteins in solution. Improved representations of solvent heterogeneity, explicit treatment of ion effects, and dynamic solvent–solute interfaces are likely to overcome limitations in modeling solvation and entropic effects.