Towards Accurate Microscopic Calculation of Solvation Entropies: Extending the Restraint Release Approach to Studies of Solvation Effects
The evaluation of the solvation entropies is a major conceptual and practical challenge. On the one hand, it is interesting to quantify the factors that are responsible for the solvation entropies in solutions, while on the other, it is essential to be able to assess the contributions of the solvation entropies to the binding free energies and related properties. In fact, the solvation entropies are neglected in almost all the studies of the binding entropies. The main problem is that widely used approaches, such as the quasiharmonic (QH) approximation do not provide reliable results particularly, in cases of shallow potential and multidimensional surfaces while brute force evaluations of the entropic effects by simulating temperature dependence of the free energy converges very slowly. This paper addresses the above issue by starting with an analysis of the factors that are responsible for the negative solvation entropy of ions, showing that it is not due to the change in the solvent vibration modes or to the solvent force constant but to the changes in the solvent configurational space upon change in the solute charges. We begin by clarifying that when one deals with aqueous solutions, it is easy to evaluate the corresponding entropic effect by the Langevin dipole(LD) treatment. However, in this work we are interested in developing a general microscopic tool that can be used to study similar effects in the proteins. To this end, we explore the ability of our restraint release (RR) approach to evaluate the solvation entropy. We start this analysis by reviewing the foundation of this approach and in particular, the requirements of minimizing the enthalpy contribution to the RR free energy. We then establish that our approach is not a specialized harmonic treatment but a rather powerful approach. Moving to the main topic of this work, we demonstrate that the RR approach provides quantitative results for the solvation entropies of monovalent and divalent ions and effectively captures the physics of these entropic effects. The success of the current approach indicates that it should be applicable to the studies of the solvation entropies in the proteins and also, in examining hydrophobic effects. Thus, we believe that the RR approach provides a powerful tool for evaluating the corresponding contributions to the binding entropies and eventually, to the binding free energies. This holds promise for extending the information theory modeling to proteins and protein-ligand complexes in aqueous solutions and consequently, facilitating computer-aided drug design.
I. Introduction<!>II. Physical Considerations<!>III. Using the RR Approach<!>III. 1. Formulation of the Restraint Release Approach<!>III. 2. Non-Harmonic Case<!>III. 3. Practical time saving by using the initial quasiharmonic treatment<!>III. 4. Optimizing the Restraint Coordinates<!>III. 5 Actual Implementation<!>IV. RESULTS<!>V. DISCUSSION<!>CONCLUSIONS AND FUTURE APPLICATIONS<!>
<p>Entropic effects play a crucial role in biological and chemical processes contributing to folding,1–6 binding processes 3,7–11 and sometimes, modulating catalytic effects.12–15 Therefore, it is important to be able to assess the magnitude of the entropic contributions in a quantitative way. Unfortunately, this task is extremely challenging and popular approaches such as the quasiharmonic (QH) and related treatments are quite inaccurate.16,17 Several approaches have been proposed in the recent years for more accurate calculation of configurational entropies1,6,12,13,17–27 and of these, some have been validated in the actual studies of biological systems.12,13,17 Furthermore, it has been shown that the entropic effects can be evaluated by brute force approaches but these require an extremely large computer time.15,17,28 Interestingly however, most of the attempts to study the entropic contributions to binding and catalysis conveniently ignored the effect of the solvation entropies.10,17,29–32 In some case s, it was found that the solvation entropy for the given process in water can be determined quite accurately by the phenomenological approach of Florian and Warshel33 who calibrated the polarization of their Langevin dipole(LD) model to reproduce the solvation entropies (see below). While the LD approach is extremely effective in solution studies, it is not clear how to implement this type of approach in the studies of the solvation entropies in the proteins. In such a case, one needs a fully microscopic approach and the development of a proper strategy which is the subject of the present work. This paper begins by considering the physical basis of the solvation entropies. It is shown that the main contributions are not associated with the solvent vibrations but with the solvent configurations. We also explore the ability of the restraint release (RR) approach to evaluate the solvation entropies. It is found that the RR approach provides an effective tool for the quantitative evaluation of the solvation energies of the charged molecules.</p><!><p>In studying the solvation entropies, it is important to have a working hypothesis regarding the origin of these effects. The simplest idea is that somehow the vibrational frequencies of the first solvation shell are higher near uncharged solute than near the solute in its charged form. One way to explore this idea is to look at the solvent vibrations that are associated with the solute-solvent interactions whilst the solute is charged and uncharged. This type of analysis can be done conveniently by the dispersed polaron (Spin Boson) method34 which provides the projection of the solvent modes along the solvation coordinates. The results of such an analysis for sodium and calcium ion in water gave the results summarized in Fig 1. Evidently, the vibrational frequencies are similar and thus, unlikely to account for the solvation entropy.</p><p>Attempts to use the calculated spectrum in the evaluation of the vibrational entropy S(v), ΔSvib =∫S (v)dv (see ref 35 for a related treatment of free energy), gave basically the same results for the charged and the uncharged systems. These results were obtained regardless of the length of the simulation or resolution of the low frequencies. Now, this finding may appear counterintuitive considering the expectations of the larger solvent force constant for both the charged and the uncharged states. Therefore, we estimated the solvent force constants by evaluating the microscopic Marcus parabolas using our well tested and widely used formulation. 36 It was found that the curvature of the charged and the uncharged parabolas are very similar (also see ref 36 for a related study). This indicated that the solvent vibrations cannot account for the solvation entropy.</p><p>The above finding tells us to look in the same direction as explored in the LD treatment. That is, it is more likely that the solvation entropies reflect the tendency of the dipoles of the first solvation shell to be oriented toward the charged solute and to assume random conformation when the solute is uncharged. This has been in fact the idea beyond our LD treatment33 which will be discussed below.</p><p>According to the LD model, the hydration entropy for charged solutes is dominated by the solvent immobilization in the presence of the solute electrostatic field, ξ. This entropy contribution can be estimated using the concept of accessible configurational volumes introduced by Frank and Evans37,38 and expressed as</p><p>
(1)ΔSimmob=∑jkBln(Vslj/Vplj),j∈solventdipoles, where kB is the Boltzmann constant while Vpl and Vsl are the accessible configurational volumes of the solvent dipoles in the pure (bulk) liquid and in the presence of the solute, respectively. The accessible configurational volumes can be expressed as the products of the volumes corresponding to the translational, rotational, and vibrational degrees of freedom of the solvent molecules. The use of the LD solvation model, in which the positions of the solvent dipoles are fixed, corresponds to the assumption that the translational contributions to the hydration entropy are negligible. Although this contribution is not negligible, it is certainly small compared to the changes in rotational degrees of freedom of the solvent molecules. Similarly, the vibrational entropy of water is supposed to be unaffected by the presence of the dissolved ions. Using the above assumptions and the fact that the direction of a solvent dipole in large electrostatic field is bound to lie between 0 and 2α (Fig. 2) one obtains the relationship</p><p>The angle α corresponds to the average angle between the direction of the solvent dip ole and the external field. In the LD solvation model, it is assumed that the average solvent polarization is determined even in solution by the Langevin function L(x)39: 
(3)cosα=L(x)=coth(x)−1/x
(4)wherex=μ0ξ3kBTC</p><p>In Eq 4, μ0, ξ, kB, C and T denote, respectively, the magnitude of the solvent dipole, the total electrostatic field at the position of this dipole, the Boltzmann constant, a scaling constant and the thermodynamic temperature. The assumption that the LD model is valid in solution40 has been shown to be physically valid 41 and thus can be used in entropy studies. At any rate, using Eqs 1–3, the immobilization entropy of the jth Langevin dipole becomes</p><p>Eq 5 does not include explicitly the dipole-dipole interaction term. Thus, this equation loses its validity for solvent dipoles positioned in the regions where the magnitude of the solute electrostatic field is comparable to the field from other solvent dipoles. Consequently, we introduced an empirical correction function</p><p>
(6)ΔScorj=−5kBarctan(x27)exp(x632) that decreases the magnitude of ΔSimmob for small fields. For +1 or −1 charged solutes, this correction affects the solvent dipoles beyond the first hydration shell, whereas for +2 and −2 charged solutes this term affects dipoles outside the first two hydration shells. Finally, a multiplicative constant, Ci = 0.878, that incorporates entropy contributions that cannot be described by the dipolar solvent model, was introduced and adjusted by comparing the calculated and experimental hydration entropies of atomic ions. The resulting empirically corrected immobilization entropy is expressed as</p><p>The behavior of the uncorrected (Eq 5) and the corrected (Eq 7) immobilization entropy as a function of x is presented in Fig. 2 of ref33. Finally, the total hydration entropies are calculated as the sum of the hydrophobic and immobilization entropies</p><p>This expression can be applied to both neutral and ionic solutes. Note that it has been suggested that the hydration entropy of large positive ions may be significantly different from that of large negative ions of the same size.42,43 In our model, this entropy does not depend on the sign of the solute charge since we believe that the magnitude of this dependence is still an open question.</p><p>It should be mentioned here that ΔSphob could have also been expressed by considering the decrease in the number of configurations available for the solvent dipoles near the nonpolar solute surface relative to the configurations available in the bulk solvent.44 A rigorous implementation of this model would require a complicated treatment of the explicit interactions between each surface molecule and its environment, which involves both the neighboring solvent molecules and the solute surface.</p><p>In order to explore the above idea, we evaluated the polarization of the water molecules around Na+, K+ and Ca 2+ions and around the uncharged forms of these ions. The corresponding results are summarized in Figs 3–6.</p><p>Figs 3 and 4 describe the fluctuations of the angle between the solvent field and the bisector of the two O-H vectors of the representative water molecules. Figs 5 and 6 provide histogram of the angle, θ. As expected, the figures indicate that the orientation of the molecules in the first solvation shell is more restricted in the charged solute than in the neutral solute. However, the trend in the case of molecules in the second and third solvation shell is less clear.</p><p>While some qualitative aspects of the origin of the entropic effect are observed in Figs 3–6, we need a more quantitative way that can also be "imported" to similar studies in the proteins. Thus, we turn our attention to the RR approach which is described in the next section.</p><!><p>The RR approach is based on the idea that with strong harmonic Cartesian restraints we have no entropy effects so, in some way, the free energy of releasing the restraints can tell us about the magnitude of the entropic effect. This general idea has been proposed by Herman and Wang45 for evaluating entropic effects in limited internal coordinates and is similar in some respects to that suggested in an early proposal.46 However, the problem is to formulate the RR idea correctly in a general set of coordinates and also, to take into account the fact recognized in ref 47, that the RR free energy can include a large enthalpic contribution. The resolution of this problem and the introduction of a practical general RR approach have been introduced in refs12,47 with a general cartesian treatment and an optimization approach that eliminates the enthalpic contribution from the RR free energy. The main points of this approach and the additional recent improvements will be discussed below.</p><!><p>In order to obtain the solute entropic contribution to a transfer between two potential surfaces (UI→UII), one has to separate the corresponding partition functions to the solute and the solvent contributions (where in the present case, the 'solute' includes explicit solvent molecules and the 'solvent' is defined as the rest of the system) by writing</p><p>
(9)QN=∫dR∫dre−UN(R,r)β=∫dR¯[∫dR∫drδ(R−R¯)e−UN(R,r)β]=∫dR¯qN(R¯)=∫dR¯e−WN(R¯)β and</p><p>
qN=∫dR∫drδ(R−R¯)e−UNβ where UN is the potential surface for state N and δ (R − R̅) indicates that the corresponding function would be collected at ±ΔR/2. Here R and r are the solute and the solvent coordinate, respectively, qN (R̅) is the solute unnormalized probability distribution evaluated at the given R̅ averaged over all the solvent coordinates and WN (R̅) is the corresponding potential of mean force (PMF) which includes, of course, the solute and solvent contributions. Expanding this potential around its minimum gives</p><p>
(10)WN(R¯)=UN(R¯)+ΔgsolN(R¯)=WN0(R¯0N)+(KWN/2(R¯)−R¯0N)2 where we could use R but instead, took R̅ as our variable to indicate that this is our restraint coordinate. Here, Δgsol (R̅) is the solvation free energy at R̅, KW is the force constant of the quadratic term of our expansion and 
R¯0N is the value of R̅ at the minimum of WN. The harmonic approximation of the PMF will be removed below.</p><p>For the sake of simplicity, we continue by examining a case where the PMF corresponds to a one-dimensional solute coordinate, although extension to many dimensions is straightforward. Our one-dimensional case is described by the thermodynamic cycle of Fig 7. In this cycle, we divide R to segments that can be defined by the delta function of Eq 9 or by introducing a strong quadratic constraint ((K′1/2)(R̅ − R̅′)2). In this case, we have</p><p>
(11)W′N(R¯)=WN0(R¯0N)+(KW/2)(R¯−R¯0N)2+(K′1/2)(R¯−R¯′N)2=WN(R¯′)+(K1/2)(R¯−R¯′N)2 where the notation W′ designates the addition of the constraint.</p><p>Next, we use Eq 11 and our quadratic constraint and evaluate the relevant partition functions.</p><p>
(12)QN=∫−∞∞e−WN(R¯)dR¯=e−WN(R¯0)β∫−∞∞dR¯e−(KWN/2)(R¯−R¯0N)2β=e−WN0(R¯0)β2π/βKWN
Q′N(R¯′)=∫−∞∞e−W′N(R¯)βdR¯=e−WN(R¯′N)β∫−∞∞dR¯′e−(K1/2)(R¯−R¯′N)2β=e−WN(R¯′N)β2π/βK1 where we use Q′ or the partition function evaluated with W′N and where 
WN0=WN(R¯0N). Using the above equation we can write</p><p>
ΔGI→II=−β−1ln(QII/QI)=ΔWI→II0−β−12ln(KWI/KWII)
(13)ΔG′I′→I(R¯′)=GI−G′I(R¯′)=−β−1ln(QI/Q′I(R¯′))=−β−12ln(K1/KWI)+(WI(R¯0I)−WI(R¯′I))≤0
E¯I=kBT2(∂lnQI/∂T)v=WI0+(kB/2)T−(T∂WI0/∂T)=WI0+(kB/2)T−TΔSsolI(R¯0I) where E̅ is the average energy and Ssol is the contribution from the surrounding of the specified system included in our explicit treatment. We also have</p><p>Now, we can write</p><p>Thus,</p><p>Finally, from Eq 13 and Eq 15 we have in the special case when 
R¯′I=R¯0I and 
R¯′II=R¯0II</p><p>The usage of this relationship however, requires knowing the values of 
R¯0I and 
R¯0II. This can be done, in principle, by minimizing the corresponding PMF (see section III. 4). Furthermore, using Eq 13 we can write</p><p>Thus, we estimate ΔG′(R0) by trying different R̅'s and taking the ΔG′(R̅) with the smallest absolute value. Now, we can use Eqs 17 and 18 to evaluate ΔS′I→II and 
ΔS′I′→II(R¯0I). The R̅ that minimizes |ΔG′| is also used in evaluating the 
ΔSsolI→II of Eq 17 and the total entropy ΔSI→II.</p><!><p>The validity of the inequality of Eq. 19 has been established above for harmonic systems but probably it is also valid for an harmonic potentials. Here, we only demonstrate the validity in some limiting cases. We first consider a square well potential (Fig. 8)</p><p>In this case, we can define a unit length ΔR by</p><p>
(20)ΔR=2π/βK1 and then we have,</p><p>
(21)QI=nΔRe−W0Iβ where n = L/ΔR and L is the width of the potential well. This leads to</p><p>Now, we also have</p><p>Using the above relationship for transition from I to II, with the same K1, we obtain</p><p>
(25)ΔSI→II≅kBln(nII/nI)=ΔGII′→II−ΔGI′→II as long as we are in the region where R ≥ 0 and W = W0, and under the condition that the contributions for ΔGI from the region with 
W=W1I are small. Once we move to the region where 
W=W1I, where 
W1I≫W0I (e.g. at R < 0), we have</p><p>Hence,</p><p>
∣ΔG′I′→I(0<R′1<L)∣<∣ΔG′I′→I(R′1<0)∣ thus, satisfying our condition of obtaining ΔSI→II from R′1 that minimizes ΔG′. The same treatment can be used for periodic systems such as the one in Fig. 9.</p><p>In this case</p><p>
(27)WI(R)=WI0+∑l(K/2)(R−Rl0)2θ(R)l where θ(R)l = 0 when R is out of the lth range, and θ(R)l = 1 inside this range. This gives</p><p>
(28)QI=∑l∫Rl−Rl−exp{−(WII0+(K/2)(R−Rl0)2)β}dR
=N∫−Δ/2+Δ/2exp{−(WI0+(K/2)(R−Rl0)2)β}dR≅Nexp{−WI0β}2π/Kβ where 
Δ=Rl+−Rl−. Now, we can write</p><p>We also, obtain</p><p>
(30)lnQ′≅WI(R′l)−β−1ln2π/βK1 where R′l is the center of the constraint potential. We also, have</p><p>Following the same consideration as in the previous case, and utilizing the fact that 
WI(R′l)>WI0 when 
R′l≠Rl0, we will obtain the same inequality.</p><!><p>The RR treatment is based on starting the simulations with a strong positional restraint where the entropic contributions are negligible. However, the free energy perturbation (FEP) release requires an enormous amount of computer time even where K1 is significantly large. Fortunately, we realized that the RR procedure is not really needed at the stage when K1 is large since the corresponding entropic contribution can be easily estimated by the QH approximations.48,49 In general, the QH tends to be valid when restraints are significant and only starts to become problematic when K1 is small (resulting in a range of very shallow and an harmonic potential energy surfaces).16,17,47 We exploit this fact by combining the QH and RR approximations by evaluating the entropic contribution at a moderate value of K1 by the QH approach followed by the RR procedure from this specific value. This treatment can be expressed by</p><p>
(32)−TΔSconf=−TΔS(K=K1)QH+min[ΔGRR(K=K1→K=0)] where the −TΔS(K = K1)QH designates the entropy computed by the QH approximation, where K1 is the initial value of the restraint and 'min' indicates the minimum value of the indicated ΔGRR. This approach has been successfully exploited in the studies of ribosome13 and more recently, a series of substituted phosphate diesters.50</p><!><p>Our next task is to optimize ΔGRR. This free energy term is obtained by a FEP/umbrella sampling approach using a mapping potential</p><p>
(33)Um=U(1)(1−λm)+U(2)λm where</p><p>
(34)U(1)=Usyst+(K(1)/2)∑j(Rj−R¯j)2
U(2)=Usyst where Usyst is the potential energy of the unconstraint system.</p><p>If we use the end point linear response approximation (LRA)instead of the full FEP treatment, we will have10</p><p>Now we have to find R0 that minimizes ΔGRR. This is done by</p><p>Now,</p><p>Thus,</p><p>We also know that the step from the current R to the optimal R is given by the Newton Rhapson relationship51</p><p>This gives</p><p>The LRA approximation of Eq 35 works the best when K(1) is small since this is when we get the largest contribution to ΔGRR. Since Eq 39 is not exact (due to the use of the LRA approach), we found it useful to scale it and to use</p><p>
(40)ΔR¯=α〈R−R¯〉U1 where α is a scaling factor.</p><!><p>All calculations were performed using the MOLARIS software package, and the ENZYMIX force field.10,52</p><p>The constraint release free energy terms are evaluated by a free energy perturbation (FEP)53,54 (also see ref55). The FEP method evaluates the free energy associated with the change of the potential surface from U(1) to U(2) by gradually changing the potential surface using the relationship</p><p>
UmÎ»m=U(1)(1âÎ»m)+U(2)Î»m where 
U(1)=Usyst+(K(1)/2)∑jRj−R¯j)2, U(2) = Usyst and λm is a mapping parameter that changes between (0 ≤ λm ≤ 1). The free energy increment, associated with the change of Um can be obtained by</p><p>
exp{−δG(λm→λm′)β}=〈exp{−(Um′−Um)β}〉m where 〈 〉m indicates that the given average is evaluated by propagating trajectories over Um and β − (kBT)−1. The total free energy difference between the states A and B is obtained by changing λm in n equal increments and evaluating the sum of the corresponding δG.</p><p>The RR calculations require imposing strong positional constraints on the Cartesian coordinates of the solvent molecules in the presence of a solute (either in its charged or uncharged state) whose position is fixed by applying a moderate positional constraint K=10 kcal mol−1Å−2 and to subsequently evaluate the free energies associated with releasing these restraints using the following equation: 
(41)−TΔSconf=min(ΔGRRchg)−min(ΔGRRunchg)</p><p>
ΔGRRchg and 
ΔGRRunchg denotes RR free energies of the systems with a fixed, charged and uncharged solute, respectively. All RR free energies contain a residual contribution from the enthalpy of the system. However, this contribution approaches zero for restraint coordinates that give the lowest RR energy.45,47</p><p>To calculate the optimized restraint coordinates (R̅) that would give the optimal ΔG for entropy calculations, the given system with the solute in its charged state, was subjected to MD simulations using a constraint force of 5.0 kcal mol−1Å−2 and a scaling factor (α) of 0.1 while averaging 〈R − R̅〉 over 1000 steps (see section III. 4). A total of 10 iterations were carried out in this procedure. Since Eq 39 is expressed in Cartesian coordinates and the application to water molecules resulted in some stretching of the bonding coordinates, we found it useful to relax the constraint coordinates obtained from Eq 40 in a 10 ps simulation run and generate nine different Rs(we report results for those with the lowest RR energy). Apparently, our exploratory studies indicated that using a small number of Rs may result in missing the one that minimizes the enthalpic contribution despite using Eq. 40. The best estimate of −TΔS is obtained by taking the corresponding values from the run with R that gives the smallest |ΔG|. This variational minimization reflects the fact that all the RR free energies contain enthalpic contributions and these contributions approach to zero for restraint coordinates that gives the lowest RR contribution.</p><p>An initial 30 ps molecular dynamics simulation was performed on the solvent molecules around the fixed charged and uncharged solute, with snapshots taken every 5 ps. MD simulations, 40 ps in length, were performed on each snapshot (with coordinates collected every 0.01 ps) to calculate the QH entropy for each snapshot, and this procedure was carried out for K= 3.0 kcalmol−1 in magnitude. This allows the calculation of the absolute entropy for this restraint. Subsequently, the free energy for the release of these restraints was evaluated using the RR approach employing Eq 32.</p><p>The RR-FEP calculations were performed with an 18Å simulation sphere of explicit water molecules subject to the surface constraint all-atom solvent (SCAAS) boundary conditions.56 The RR simulations, consisting of 42 windows each with a simulation time of 40ps, involved the release of the position restraints in 5 FEP stages, changing K1 from 3.0 to 1.0, from 1 to 0.3, from 0.3 to 0.03, from 0.03 to 0.003 and finally, from 0.003 to 0.0003, where all the values of K1 are given in kcal mol−1 Å−2. These simulations were done at 300K with 1fs time steps.</p><p>Simultaneously, the LD model33 implemented in the program ChemSol 2.057 was also used to evaluate the solvation entropy. This method evaluates the solvation entropy by estimating the restriction of the solvent motion due to the field from the solute.</p><!><p>As stated above, we tried to explore in this work the ability of the RR approach to reproduce solvation entropies. The most obvious benchmark is provided by the solvation entropies of ions and this challenge has been addressed here by evaluating the solvation entropies of Na+, K+ and Ca2+ ions. The corresponding RR calculations were performed on the inner 29 water molecules that are closest to the ion, while the rest of the 18Å SCAAS system to move with a very small position constraint on the oxygen atoms in addition to the SCAAS polarization and distance restraints. In short, the entropic effect associated with the effect of the solute on the inner 29 water molecules was explored while representing rest of the system by the SCAAS model. In calculating the solvation entropy, one should also consider the hydrophobic energy of forming the cavity for the uncharged ion. This hydrophobic contribution is small for the ions considered here and were added to the final calculated results. This was done using the corresponding values obtained from ChemSol. In future, we hope to obtain such contributions from our RR calculations. The results of these calculations are summarized in Tables 1–3. Remarkably, we reproduce the observed trend in an almost quantitative way. The tables also give the quantitative LD results but as mentioned before, the LD treatment is empirical in nature and cannot be used within the protein active site. It should also be mentioned that the recent works of Aqvist and coworkers28 reproduced quantitative results for the same system studied here. However, this was done by evaluating the van't Hoff approach that requires a major investment in terms of computer time.</p><p>The fact that the computational methodology employed in this work reproduces the observed solvation entropy indicates that the physics of this effect is captured by this approach and that the effect is associated with the change in the effective space available for the solvent molecules upon change in the solute charges.</p><p>The overall trend is clearly discussed in section II, but here we have a well-defined simulation approach that reproduces the entropic effect quantitatively. Furthermore, in order to clarify the nature of the challenge overcome by the RR calculations, we also report in Table 4 the QH results obtained with no restraint (K= 0 kcal mol−1 Å−2) for Na+ and Ca2+ ions. Although QH entropic estimation for Na+ ion seems fairly accurate, it provides completely erroneous results for Ca2+ ion. Thus, we may say that the QH is simply unable to capture the solvation entropies. This is probably one of the best illustrations of the fundamental problems with the widely used QH and related approaches.</p><!><p>The present work focuses on the calculations of the entropic effects, despite the fact that most questions about biophysics functions are related to the corresponding free energy surfaces.58,59 Additionally, we are well aware of the fact that free energy calculations converge much faster than the entropic calculations. However, in many cases complete understanding of the entropy-enthalpy compensation is useful for the deeper understanding of the underlying physics of the system and in particular, the origin of the temperature dependence of various processes. Moreover, the evaluation of the entropic effects is one of the utmost challenges of computational biophysics and in fact, the failure to reproduce entropic effects is frequently used by the experimental community in trying to diminish the important free energy calculations (also see discussion in ref60).</p><p>In any case, the accurate evaluation of the entropic contributions to the free energy changes in condensed phases is an extremely challenging task. Popular approaches such as the QH and related treatments have been widely used to extract the configurational entropy from molecular dynamics simulations.48,61,62 Unfortunately, although QH approximation functions quite well for simple systems with single energy wells, this approach provides a poor approximation of the configurational entropies of systems with multiple occupied energy wells, as it merges multiple narrow energy wells into one broad energy well.16,17 In fact, the validity of such approaches in systems with very shallow potentials is questionable. Therefore, the use of the QH approach in the studies of binding energies or related problems is problematic. Several approaches have been proposed for obtaining quantitative entropic contributions, but most of these approaches have not been able to provide quantitative results in general cases, except by brute force studies of the free energy change with temperature that requires enormous simulation time for convergence.28 In our view and experience, one of the most promising strategies is our RR approach that has been shown to give quantitative results for activation entropies.10,12,13,47,50 Nevertheless, this approach has not been validated in the studies of binding entropies.</p><p>One of the most serious problems in calculations of binding entropy is the evaluation of the solvation entropies in water and in the protein active site. This contribution has been generally neglected in key studies29–32 although it has been evaluated or included in calculations of activation free energies.13,15,47 The problem of activation entropies is also a problem of general interest both for studies of molecules in solutions and in proteins.</p><p>When one deals with molecules in solution, it is possible to use the LD models to obtain a relatively reliable estimate of the solvation entropies. Unfortunately, this approach has not been generalized to protein active sites where one must use a microscopic approach. The studies of ions in solution should provide the best benchmark for the development and validation of the corresponding approach since the corresponding values are known experimentally.</p><!><p>The molecular simulations of the entropy contribution in the areas of protein-ligand binding and enzyme catalysis, is an extremely challenging task and the subject of intense research. The analyses of the role of solvent entropies in ligand binding and enzyme catalysis are still rare, although knowledge of their values would be of great utility considering that the binding processes can be relatively dominated by solvent entropy effects. Therefore, there is a strong motivation to obtain the solvation entropic contributions that can lead to a better understanding of the process under consideration. Although the solvent entropies can be extracted from the temperature dependence of the explicit solvent FEP simulations, these calculations are difficult to converge and their convergence tends to be considerably slower than that of the free energy itself. The present work focuses on the development of efficient microscopic model for the studies of the solvation entropies in aqueous solution. Our preliminary studies clarified that the solvation entropies are not associated with the solvent vibrations but rather with the availability of larger configurational space upon transition from the charged to the uncharged form of the solute. This effect is easily observed in the cases of divalent ions but in the case of monovalent ions it is restricted to the first solvation shell. With this in mind, we demonstrated that the solvent entropies can be successfully and accurately evaluated by our RR approach. This has been illustrated for the simple but relevant test case of monovalent and divalent metal ions. It is clear from the results that solvation entropies depend strongly on the magnitude of the charge of the solvated ion. This is obviously related to the corresponding restriction in the solvent fluctuations.</p><p>Although convergence problems may exist for larger and more complex simulation systems, a good agreement is obtained between the calculated and the experimental data for solvation entropies of monovalent and divalent metal ions. Furthermore, we also developed a systematic approach for the minimization of ΔGRR through the Newton-Rhapson minimization of the ΔGRR free energy, formulated within the LRA. Our previous studies have shown that the random search for the optimal R̅ may otherwise, require extensive simulations. Thus, this approach reduces computer cost and time while finding the optimal R̅ for our RR calculations. Overall, the use of the RR approach has been found to provide a powerful way of quantifying the solvation entropies. On the other hand, it has been demonstrated that the QH approach provides an entirely inadequate tool for this important purpose.</p><p>The success of the present approach in the evaluation of the electrostatic components of the solvation entropies indicates that the RR method will probably be useful in the studies of hydrophobic effects where the physics (of restricted orientation) is smaller. Efforts are underway in our laboratory to extend the utility of this approach in the studies of the contribution of the solvation entropies to the binding entropies, hoping to help advance the understanding of this important effect. The corresponding strategy will be similar to the one used in the present study and will require releasing the restraints on the protein residues and the water molecules within a given cut-off distance. However, this approach will require careful refinement since releasing the restraints on all of the protein residues may lead to computational instability while releasing the restraints on a few residues may not capture the full solvation entropy contribution. In future, it should be feasible to estimate the magnitude of the solvent entropic contributions involved in ligand binding to proteins. Generally, the binding affinity of a compound can be improved by generating a favorable binding enthalpy, favorable solvation entropy, or by minimizing the unfavorable conformational entropy. In some cases, it is much easier to optimize the entropy. For example, medicinal chemists have attempted to conformationally constrain the substrate thus, restricting the entropy loss upon binding. Further improvement may be obtained by trying to optimize the solvation entropies by choosing specialized functional groups. However, the main advancement in our approach is in obtaining more reliable predictions of the entropic contributions to ligand binding. Our RR approach should therefore, be viewed as a step towards the generation of algorithms capable of thermodynamic dissection particularly, those that can predict the entropic consequences of introducing different chemical functionalities in the molecular scaffold generated either de novo (during lead generation) or while being optimized (during lead optimization)for improving their binding affinity towards a target protein. Thus, this study may have implications towards our understanding of protein-ligand interactions, for theoretical modeling of molecular recognition, and for rational drug design.</p><!><p>The Fourier transform of the time dependence of the energy gap (ΔVab) between the charged (state a) and uncharged (state b) states of (a) Na+ ion (b) Ca2+ ions. The dimensionless coordinate displacement (Δ) values are obtained from a simple Fourier transform on the time dependent energy gap obtained from classical trajectories on the potential surface of Va.</p><p>Schematic representation of the restricted configurational volume of a grid-centered solvent dipole in the presence of an external electrostatic field, ξ. In the pure liquid, the rotation of the solvent dipole is assumed to be free and the accessible volume corresponds to the outer shell of a sphere of a radius μ0. In the presence of the external field, the end point of the solvent dipole becomes confined, on the average, to the outer shell determined by the solid angle between the direction of the field and the maximum deviation of the solvent dipole from this direction.</p><p>The time-dependent fluctuations of three water molecules at distances of 1.96Å (a), 4.90Å (b) and 6.21Å, respectively (c) from charged (left) and uncharged (right) Na+ ion during a 2ns molecular dynamics simulation.</p><p>The time-dependent fluctuations of three water molecules at distances of 1.96Å (a), 4.90Å (b) and 6.21Å, respectively (c) from charged (left) and uncharged (right) Ca2+ ion during a 2ns molecular dynamics simulation.</p><p>Histogram of the average orientation of three water molecules at distances of 1.96Å (a), 4.90Å (b) and 6.21Å, respectively (c) from charged (left) and uncharged (right) Na+ ion during a 2ns molecular dynamics simulation.</p><p>Histogram of the average orientation of three water molecules at distances of 1.96Å (a), 4.90Å (b) and 6.21Å, respectively (c) from charged (left) and uncharged (right) Ca2+ ion during a 2ns molecular dynamics simulation.</p><p>A schematic illustration of the effect of R̅ and the enthalpic contribution to ΔG′. The figure considers one-dimensional potentials of mean force (W) for two systems (e.g. charged and uncharged solute) and the effect of confining the system to small segments of R by a strong constraint (K1). The figure also represents the effect of restraining the system in different R̅'s by the PMF of the restraint systems, (W′). The value of ΔG′ depends on R̅ and the ΔG′ obtained with the minimum of the corresponding W (i.e., R̅0) has the smallest absolute value. This ΔG′ is our −T ΔS′.</p><p>Schematic representation of a square well potential.</p><p>Schematic representation of a periodic system that may be utilized to obtain ΔSI→II from R′1 that minimizes ΔG′.</p><p>Entropy calculation for a solvated Na+ ion, using the RR approacha</p><p>Energy values are given in kcal/mol. The table includes the results of six sets (1–6) of simulations with different restraint coordinates. −T(ΔS)hyd is the cavity formation contribution calculated using Chemsol. The value in bold signifies the smallest ΔG and thus, the value that satisfies Eq 32.</p><p>Entropy calculation for a solvated K+ ion, using RR approacha</p><p>Energy values are given in kcal/mol. The table includes the results of six sets (1–6) of simulations with different restraint coordinates. −T(ΔS)hyd is the cavity formation contribution calculated using Chemsol. The value in bold signifies the smallest ΔG and thus, the value that satisfies Eq 32.</p><p>Entropy calculation for solvated Ca2+ ion, using RR approacha</p><p>Energy values are given in kcal/mol. The table includes the results of six sets (1–6) of simulations with different restraint coordinates. −T(ΔS)hyd is the cavity formation contribution calculated using Chemsol. The value in bold signifies the smallest ΔG and thus, the value that satisfies Eq 32.</p><p>QH calculations for solvation entropiesof the metal ions.a</p><p>Energy values are given in kcal/mol. The table includes the results of six sets (1–6) of simulations with K=0 kcal mol−1 Å −2 on the system. The final entropy value is obtained by taking an average of the difference between the charged and uncharged form of the six different conformations used.</p>