Why Solvent Response Contributions to Solvation Free Energies Are Compatible with Ben-Naim’s Theorem

We resolve a seeming paradox arising from a common misinterpretation of Ben-Naim’s theorem, which rests on the decomposition of the Hamiltonian of a molecular solute/solvent system into solute–solvent and solvent–solvent interactions. According to this theorem, the solvation entropy can also be decomposed into a solute–solvent term and a remaining solvent–solvent term that is commonly referred to as the solvent reorganization term. Crucially, the latter equals the average solvent–solvent interaction energy such that these two solvent–solvent terms cancel and thus do not change the total solvation free energy. This analytical result implies that changes in the solvent–solvent interactions cannot contribute to any thermodynamic driving force. The solvent reorganization term is often identified with the contribution of many-body solvent correlations to the solvation entropy, which seems to imply that these correlations, too, cannot contribute to solvation. However, recent calculations based on atomistic simulations of a solvated globular protein and spatially resolved mutual information expansions revealed substantial contributions of many-body solvent correlations to the solvation free energy, which are not canceled by the enthalpy change of the solvent. Here, we resolved this seeming contradiction and illustrate by two examples—a simple Ising model and a solvated Lennard-Jones particle—that the solvent reorganization entropy and the actual entropy contribution arising from many-body solvent correlations differ both conceptually and numerically. Whereas the solvent reorganization entropy in fact arises from both solvent–solvent as well as solute–solvent interactions and thus has no straightforward intuitive interpretation, the mutual information expansion permits an interpretation in terms of the entropy contribution of solvent–solvent correlations to the solvation free energy.


Introduction
The hydrophobic effect is an essential driving force for many processes in nature, such as phase separation, membrane formation, [1][2][3] or the function and folding of proteins. 4,5Despite its significance, the hydrophobic effect is not yet fully understood from first principles, and hence its molecular explanation remains controversial. 6The early "iceberg hypothesis" by Frank and Evans, 7 for example, turned out to be equally popular and controversial.Frank and Evans explained the unfavorable solvation free energy of hydrophobic solutes in water by an entropic penalty due to an ordered "iceberg" structure of water molecules that forms around the solute.The term "iceberg" is not meant to be taken literally, but rather refers to a higher ordering of the first few solvation shells compared to bulk water. 7,8ndeed, the hydrophobic effect has been ][10][11] However, in seminal papers Ben-Naim [12][13][14] and Yu et al. 15 have analytically proven that for pairwise interactions, the entropic and the enthalpic parts of the waterwater interactions exactly cancel and, therefore, do not contribute to the net free energy change of solvation.This important theorem has led to the general understanding that any solventresponse to the solute, such as the Frank and Evans "icebergs", has a net-zero effect on the free energy change and that, therefore, the solvent response cannot drive solvation.However, ordered structures of water molecules around hydrophobic solutes are still implied as the cause of hydrophobicity, 16,17 a notion that has then been corrected by others. 14,18Overall, this seeming contradiction caused considerable confusion, and still does.Against this background, our recent finding that water-water correlations contribute markedly to the folding free energy of crambin 19 bears the potential for further confusion, as indeed testified by anonymous referee reports.Their comments have prompted us to look deeper into the subject, and to share our analysis as well as two illustrative examples with a broader readership.
Specifically, in our previous paper 19 we calculated and compared solvation shell entropies of the solvated folded conformation of crambin and a molten-globule-like conformation of this prototypic globular protein.To this end the method Per|Mut 20,21 was used, which employs a mutual information expansion [22][23][24][25] into single-molecule entropies ∆S 1 and entropy contributions ∆S ≥2 arising from correlations between pairs and triples of -mostly nearby -solvent molecules.In our molecular dynamics simulations, the molten-globule-like conformation of crambin showed many hydrophobic residues, which are buried within the folded conformation, exposed to the solvent.Relative to the native fold, we observed indeed a marked entropic free energy contribution to the solva-tion free energy due to strongly correlated water molecules in the innermost solvation shells.This result may seem to be incompatible with the Ben-Naim theorem, in which case the notion that the solvent response cannot be a thermodynamic driving force would need to be reconsidered.
Here we will show that this finding is in fact -counter-intuitively -perfectly compatible with the Ben-Naim theorem.Our analysis will, further, provide a deeper understanding of the contribution of the solvent response to the solvation free energy.We will illustrate our reasoning by a simple Ising model example that can be exhaustively enumerated as well as by a more realistic example of a Lennard-Jones particle solvated in liquid argon.

Canonical decomposition
Ben-Naim [12][13][14] has proven that the change of average solvent-solvent interaction energies upon solvation is exactly compensated by a corresponding entropy change, such that there is no net free energy contribution.Later, Yu et al. 15 obtained essentially the same result by considering a solvation process described by the coupling parameter λ (λ = 0: not solvated, λ = 1: fully solvated).
In particular, they demonstrated that for a Hamiltonian consisting of pairwise solute-solvent (uv) and solvent-solvent (vv) interactions, the internal energy (∆U ) and entropy (∆S) changes can be expressed as Whereas the internal energy and entropy parts ∆U uv and T ∆S uv only contain solute-solvent interactions H uv , the remaining terms are referred to as "solvent-solvent" terms, ∆U vv and T ∆S vv , respectively.The important finding by Ben-Naim and Yu et al. is that these two terms are identical and thus cancel in the net free energy difference which thus only contains the so called solutesolvent terms.Note, however, that the ensemble -and thus all averages, including the solvent-solvent terms -are affected by both solvent-solvent as well as solute-solvent contributions.Due to this very fact, the terminology "solute-solvent" vs. "solvent-solvent" is highly misleading and, as will become clear further below, is the root of longstanding and widespread confusion.
In particular, the seeming absence of solventsolvent terms in the free-energy balance has led to the widely held belief that the solvent response to the presence of a solute, e.g., solvent rearrangements such as the Frank and Evans "icebergs", cannot contribute as a thermodynamic driving force. 14,26

Mutual information expansion
Additionally, because solvent-solvent terms affect the ensemble itself and therefore also -implicitly -contribute to T ∆S uv , it is also misleading to split entropic contributions, which are inherently non-pairwise, into contributions that are defined via pairwise interaction energies, as in equation 2. As an alternative, and to gain physical insight into the solvation process that can be interpreted in a more straightforward manner, we suggest to use a mutual information expansion (MIE), [22][23][24][25] as, e.g., used in the recently developed method Per|Mut. 20,21ccordingly, the total solvent entropy is decomposed into single-body entropies, akin to an ideal-gas term, and multi-body correlations, where S 1 (i) are the single-body entropies of the (three-dimensional) probability distributions of molecules 1, . . ., N ; I 2 (j, k) and I 3 (l, m, n) are the two-body and three-body mutual information terms of molecule pairs and triples, respectively.These terms are defined as and represent the entropy change due to two and three-body correlations, respectively.In this notation, S(j, k) and S(l, m, n) are the entropies of the (six and nine-dimensional) marginal distributions of the full configuration space density ϱ with respect to molecule pairs j, k and triples l, m, n, respectively, e.g., A full MIE up to the N -body correlation term yields an exact entropy decomposition.In our numerical approach below, evaluation of the respective integrals would require sampling over the full 3N -dimensional configuration space, however, which is impractical.We therefore truncated the expansion after the three-body correlations to obtain a good approximation of entropy, neglecting higher-order terms.For short-ranged interactions, these have indeed been demonstrated to be small. 27Methods

Ising model
To assess the solvent response to a solute (e.g., a protein), as sketched in Fig. 1A, we first considered the simple 4 × 4 sub-critical Ising model sketched in Fig. 1B.In this model, each spin σ i,j = −1, +1 interacts with its nearest neighbors with an interaction strength J = 0.2 under periodic boundary conditions.Here, the spins mimic a solvent with the four most center spins (shaded in red) interacting with an external field λ, which mimics the interaction with a solute.Note that, because the solute is described purely by these interactions, the Ising model depicted in Fig. 1B does not contain any explicit solute degrees of freedom.
Accordingly, the Hamiltonian reads where the first sum runs over all nearest neighbors, and the second sum runs over the spins shaded in red (i.e., the "solvation shell", the solute is not shown in Fig. 1B).The probability of each state x ∈ X = [−1, +1] 4×4 reads with the partition function Z chosen such that The entropy S and the average solvent-solvent interaction energy therefore read Following Yu et al., 15 all solute-solvent and solvent-solvent entropy changes were calculated according to equation 4.These values are therefore subject to a small integration error due to the required numerical integration, for which we used 251 discrete λ-intermediates.For the Ising model, all calculations were carried out with unitless energies, i.e., k B = T = β = 1.

MD simulations
To quantify the response of an argon-type liquid to a Lennard-Jones solute, two systems were simulated, an unsolvated system containing 512 argon-type atoms and a solvated system with an additional immobilized Van-der-Waals sphere as a 'solute'.All molecular dynamics (MD) simulations were carried out using the software package Gromacs 2020.6 [28][29][30][31][32] with a leapfrog integrator with a 2 fs time step.][35] The Lennard-Jones 36 parameters (particle size σ and potential depth ϵ) for the model solute were chosen as twice as those of argon to enhance the statistical significance of average energy differences.All Van-der-Waals interactions were switched between 1.0 nm and 1.2 nm and no dispersion correction was applied.To immobilize the solute at the center of the simulation box, the freeze-options within Gromacs were used.During all simulation runs, the temperature was kept at 120 K using the V-rescale thermostat 37 with a time constant of 0.1 ps.The unsolvated (pure argon) system was equilibrated at 1 bar pressure in a 20 ns NPT-run using the Berendsen barostat, 38 resulting in a (3.073 nm) 3 cubic simulation box.The second, solvated system was prepared by adding the solute to the simulation box and allowing for a further equilibration, lasting 10 ns under NVT conditions.For both systems, production runs, each lasting 4 µs, were carried out under NVT conditions.For subsequent analysis, configurations were stored every 10 ps, resulting in trajectories consisting of 4 • 10 5 frames each.

Entropy calculation
Entropy contributions were calculated using the method Per|Mut, 20,21 which utilizes a permutation reduction 39,40 and a mutual information expansion (MIE) [22][23][24][25] into one-, two-and three-body correlations.For permutation reduction, 50 different simulation snapshots were randomly selected as reference structures and a MIE was carried out using each of the permutationally reduced trajectories.In the MIE, the mutual information between all pairs of argon atoms was taken into account; triple-wise mutual information terms were cut off at an average distance of 0.5 nm after permutation reduction. 20,21All MIE orders were calculated using a k-nearest-neighbor algorithm with a value of k = 1.
From the resulting entropy difference ∆S MIE between the unsolvated system and the solvated system, the free energy difference ∆F MIE = ∆U − T ∆S MIE was calculated, where the internal energy difference ∆U was obtained directly from the average interaction energies in the simulation runs.
Solute-solvent and solvent-solvent entropy differences ∆S uv and ∆S vv , respectively, were calculated by thermodynamic integration (TI) from the unsolvated state to the solvated state using 200 equidistant windows, each lasting 200 ns.As a control for the Per|Mut results, the total entropy difference ∆S TI = ∆S uv + ∆S vv and the free energy change were calculated using standard TI.Errors of the internal energies were calculated as σ U / N f − 1, where σ U are the standard deviations of the respective interaction energies from N f = 400 • 10 3 simulation frames.Due to the long interval of 10 ps between frames, these were considered statistically independent.Similarly, Per|Mut errors were estimated as the standard errors resulting from the 50 permutationally reduced simulation trajectories.TI errors were estimated from the difference between two independent sets of TI simulation runs with identical input parameters but different initial (random) velocities, but turned out to be negligible for all further analyses.

Results and discussion
To investigate the seeming contradiction between the Ben-Naim theorem and the free energy effects of increased solvent correlations observed for the crambin 19 solvent shell, we calculated the relevant contributions of the solvent response to the solvation free energy for two simple model systems, for which sampling errors can be neglected.Specifically, we will compare all free energy contributions of the Yu et al. decomposition (equations 3-5) with the mutual information expansion (equation 1).
We will first consider an idealized solvation process for an Ising model, for which all relevant quantities can be exhaustively enumerated, such that the results are exact to numerical precision.Subsequently, we will consider a liquid argon-type Lennard-Jones system with a van-der-Waals solute, the relaxation times of which are short with respect to simulation times, such that for this more realistic system sampling errors can be assumed to be very small with respect to the relevant energy and entropy differences.

Ising model
As a simple illustrative model of a solute in a solvent (Fig. 1A), we consider the 4 × 4 subcritical Ising model shown in Fig. 1B.Here, each spin represents a solvent molecule that interacts with its nearest neighbors.The effects of a solute are modeled by an external field with strength λ that acts on the "solvation shell" (red), consisting of the four spins at the center.
Figure 2 shows the exact relevant thermodynamic quantities as a function of the coupling parameter λ, calculated by full enumeration as described in section 3.1.As expected, with increasing coupling to the solvent, the total entropic free energy contribution −T S TI (dotted black line) becomes less favorable (i.e., it increases), and eventually saturates at an entropy difference of −T ∆S TI = 2.48 between fully sol-vated (λ = 1) and fully decoupled (λ = 0).This contribution is dominated by the unfavorable solute-solvent contribution −T ∆S uv = 3.70 (green solid line), which is partially compensated by the favorable solvent-solvent contribution −T ∆S vv = −1.22(green dashed-dotted line).
As shown by the solid red line in Fig. 2, the average solvent-solvent interaction energies U vv increase for increasing coupling parameter λ and thus contribute unfavorably to the free energy change by ∆U vv = 1.22.Fully in line with the Ben-Naim theorem (equation 5), U vv is indeed precisely compensated by −T S vv , such that ∆U vv − T ∆S vv = 0 and, hence, ∆F = ∆U uv − T ∆S uv (∆U uv = −15.98,not shown in Fig. 2).
Does this finding imply that solvent-solvent correlations do not contribute to the solvation free energy?To answer this question, consider the above mutual information expansion of entropy, which directly quantifies these correlations.The single-body term −T S 1 (light blue line) underestimates the entropy on average by 0.53 units and contributes −T ∆S 1 = 2.82 to the overall entropy change.Inclusion of the two-and three body correlation terms (−T ∆S MIE , solid blue line) improves the approximation markedly, with an average deviation from the exact values below 0.13 units.The correlation terms −T S ≥2 (dashed-dotted light blue line) add a small favorable contribution of −T ∆S ≥2 = −0.07 to the overall MIE entropic free energy change of −T ∆S MIE = 2.75.
Crucially, the term −T S vv differs from the solvent-solvent correlations −T S ≥2 both by definition and, indeed, also numerically as shown in Fig 2 .As a result, also −T S uv differs from −T S 1 , and the entropy change due to solvent correlations −T ∆S ≥2 is also not compensated by any canonical internal energy term.This simple example illustrates that, generally, solvent correlations do contribute to the solvation free energy; it also clarifies why this finding is not in conflict with the Ben-Naim theorem.Here, spins (gray arrows) represent water molecules.The four spins shaded in red interact with the solute; this interaction is described by an external field λ that acts on the four spins.(C) As a more realistic model system, argon-type atoms (red) are coupled to a solute Lennard-Jones sphere (yellow).On the left (λ = 0), the solute is decoupled from the argon solvent; on the right (λ = 1), the solute is fully solvated.

Argon
Is this subtle but important distinction between −T S vv and the actual many-body contribution to the solvation entropy, −T S ≥2 , also relevant for more realistic systems?To address this question, we carried out MD simulations of a system comprising 512 argon-type atoms and an immobilized Lennard-Jones "solute", as described in section 3.2 (see also Fig. 1C).Here we calculated the free energy change of solvation, as well as the relevant enthalpic and entropic contributions using both, Per|Mut and thermodynamic integration.
As shown in Fig. 3, the internal energy change ∆U upon solvation is favorable and totals −9.2 kJ•mol −1 , to which solvent-solute interactions (∆U uv ) contribute −9.5 kJ•mol −1 and solvent-solvent interactions (∆U vv ) contribute In line with the Ben-Naim theorem, the later contribution is exactly compensated by −T ∆S vv = −0.3kJ•mol −1 , which, also for this system, might suggest that the solventsolvent interactions and correlations, taken together, do not contribute to the solvation free energy.
To test our assumption that four-body and higher correlations not included within −T ∆S MIE are sufficiently small, we have also calculated the relevant entropy terms using TI (green).Indeed, the similar total entropy change of −T ∆S TI = 6.6 kJ•mol −1 supports this assumption and shows that the contribution of the higher correlations to the solvation entropy is markedly smaller than the MIE estimate.Also for this more realistic system, the entropy change due to solute-solvent interactions (−T ∆S uv = 6.9 kJ•mol −1 ) dominates, and −T ∆S vv does not even describe the correct sign of the actual solvent-solvent correlation contribution to the solvation free energy.The remaining difference of ca. 1 kJ•mol −1 between the MIE and TI solvation entropies is also reflected in the respective total free energies ∆F MIE = ∆U − T ∆S MIE = −1.5 kJ•mol −1 and ∆F TI = −2.5 kJ•mol −1 , respectively, underscoring that this difference is mainly due to the truncated MIE expansion rather than sampling uncertainties.
Similar to our findings for the above Ising model, also for the more realisic argon-type system the two possible entropy decompositions differ significantly.Whereas the small size of the two solvent-solvent terms ∆U vv and −T ∆S vv -and in particular their mutual cancellation -seem to show that the solvation of this Lennard-Jones particle is unaffected by the reaction of the solvent, the actual solventsolvent entropy contributions are substantial and not compensated by any canonical internal energy term.We conclude that also for the solvation of a Lennard-Jones particle in a Lennard-Jones fluid, the induced solvent reorganization contributes markedly to the solvation free energy.

Conclusions
We pointed out that the entropy decomposition by Ben-Naim and Yu et al. into a contribution S uv from solute-solvent interactions and a remaining contribution (S vv ) as defined in equation 4, differs conceptually from direct evaluation -e.g., via a mutual information expansion -of solvent-solvent correlation contributions to the solvation free energy.In particular, the term "solvent-reorganization entropy" for ∆S vv is highly misleading, because it creates the wrong impression that any solvent response to the presence of a solute cannot contribute to the net free energy.
Two examples served to illustrate the solu-tion of this seeming paradox.First, a simple semi-analytical Ising model, which permitted exhaustive enumeration, establishes that the conceptual difference between ∆S vv and ∆S ≥2 actually gives rise to marked numerical differences.Second, our MD simulations of solvation within a Lennard-Jones liquid show that this distinction is also relevant for a more realistic solvation system.For both systems, ∆S vv is exactly compensated by the change of average solvent-solvent interactions (∆U vv ), as required by Ben-Naim's theorem.
In more general terms -as already pointed out by Lee 41 -one can define for any change of an entropy component T ∆S xy an appropriate ∆U xy component such that T ∆S xy = ∆U xy .However, such construction does not necessarily allow for a physically meaningful interpretation; in particular, it does not support the conclusion that 'xy' is irrelevant for the solvation process.Whereas the canonical decomposition of pairwise interaction energies into solvent-solute (H uv ) and solvent-solvent (H vv ) terms, as well as the corresponding decomposition of internal energies, are certainly physically meaningful, this does not necessarily apply to the corresponding entropy terms due to their inherently non-pairwise nature.Instead, an entropy decomposition into a single-body term and multibody correlations provides a more intuitive understanding.
We hope our explanations and examples will contribute to resolving a long-standing controversy and the resulting widespread confusion.Fully in line with Ben-Naim's theorem, solventsolvent correlations can -and generally docontribute markedly to the overall free energy of solvation, thus underscoring the need for an improved understanding of the "iceberg"-type ordering of solvent shells, in particular near complex macromolecular solutes and surfaces.

Figure 1 :
Figure 1: (A) Sketch of a solute (black) in water (gray angles), where the solvation shell (red area) interacts directly with the protein.(B) Sketch of a 4 × 4 Ising model, serving as a simplified model of solute-solvent / solvent-solvent interactions.Here, spins (gray arrows) represent water molecules.The four spins shaded in red interact with the solute; this interaction is described by an external field λ that acts on the four spins.(C) As a more realistic model system, argon-type atoms (red) are coupled to a solute Lennard-Jones sphere (yellow).On the left (λ = 0), the solute is decoupled from the argon solvent; on the right (λ = 1), the solute is fully solvated.

Figure 2 :
Figure 2: Thermodynamic quantities of the Ising model as a function of the external field (λ).The precise entropy contribution (−T S TI ) is shown as a dotted black line.The singlebody entropy (−T S 1 ) from the Per|Mut MIE is shown as a solid light blue line; the contribution from two-and three-body correlations (−T S ≥2 ) is shown as a dashed-dotted light blue line.The entropy from the MIE approximation, including correlations, is shown in dark blue.The entropy decomposition into −T S uv (solid line), and −T S vv (dashed-dotted line) is shown in green.The red line represents the solvent-solvent interaction energies (U vv ).For a better visual representation, U vv and −T S ≥2 are shifted by 6.61 and 8.71 units, respectively.

Figure 3 :
Figure 3: Solvation free energy contributions (in kJ•mol −1 ) of the fixed Lennard-Jones solute in an argon-type liquid.Red bars denote the internal energy change and its contributions.Blue and green bars denote the entropy change and its contributions, calculated using Per|Mut and TI, respectively.Purple bars show the overall free energy change, as calculated using Per|Mut and TI.Estimated sampling uncertainties are shown as small black bars.