System size dependence of hydration shell occupancy

TL;DR

This paper investigates how the size of the system influences the free energies associated with hydration shells, revealing that both hydrophilic and hydrophobic effects vary with system size, impacting macromolecular hydration understanding.

Contribution

It demonstrates the system size dependence of hydration free energies and clarifies the balance of hydrophilic and hydrophobic contributions within the quasichemical framework.

Findings

Hydration shell evacuation free energies depend on system size.

Hydrophilic and hydrophobic contributions decrease as system size increases.

Large systems are necessary for accurate macromolecular hydration modeling.

Abstract

The free energies to evacuate the first hydration shell around a solute and a cavity defined by the first hydration shell depend on the system size. This observation interpreted within the quasichemical theory shows that both the hydrophilic and the hydrophobic contributions to hydration depend on the system size, decreasing with increasing system size. Although the net hydration free energy benefits somewhat from the balancing of hydrophilic and hydrophobic contributions, a large system still appears necessary to describe the effect of the solvent on the macromolecule.

Equations

$${\mu }^{(\mathrm{ex})}=k_{\mathrm{B}}T\ln x_0-k_{\mathrm{B}}T\ln p_0+{\mu }^{(\mathrm{ex})}(n=0).$$

Full text

Free energy to evacuate the inner-hydration shell of volume $v_{o}$ in a cubic simulation cell of volume $L^{3}$.

How many water molecules are needed to simulate a hydrated protein? On the basis of their studies on hemoglobin, Karplus, Meuwly, and coworkers1 have found that a large system size than conventionally used — about 24 water molecules per protein heavy atom, or about 105,000 water molecules in total — is required to capture the dynamics of deoxy-hemoglobin (PDB: 2DN2). While the statistical resolution of the results and the suitability of the analysis to address system size effects have recently been questioned2, the observation in Ref. 1 suggested to us a possible reason why our free energy calculations in conformational switching of protein $G_{B}$3 proved inconclusive. In pursuing this line of research, we have uncovered a hitherto unexpected feature of finite size effect that we discuss below.

In molecular simulations of hydration, using a finite number, $N$, of solvent molecules is unavoidable. In periodic simulations, there are implicit system size effects due to periodicity, for example in pair-correlations 4 and on electrostatic self-interaction 5, 6. The are also explicit finite size effects that arise from an ensemble dependence 7, 8, for example, in calculating the compressibility 9, 10, 11. From the vantage of quasichemical theory theory 12, 13, 14, 15, we reasoned that such explicit finite size effects should be present in simulations of hydration as well.

To anchor the discussion, first consider the hydration of a simple solute, an imidazole ring (Fig. 1).

Once we demarcate the inner-shell domain, the hydration free energy of the solute is given by quasichemical theory 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22 as

[TABLE]

In Eq. 1, $-k_{\mathrm{B}}T\ln x_{0}$ is the work done to evacuate the inner shell in the presence of the solute and $-k_{\mathrm{B}}T\ln p_{0}$ is the corresponding quantity in the absence of the solute; $x_{0}(p_{0})$ is the probability to observe an empty inner shell in the presence (absence) of the solute. $k_{\mathrm{B}}T\ln x_{0}$ is termed the chemistry contribution, since it reflects the contribution to the hydration from short-ranged solute-solvent interactions. It is a measure of the hydrophilic contributions to hydration. $-k_{\mathrm{B}}T\ln p_{0}$ is the packing contribution and is a measure of primitive hydrophobic effects23, 24. $\mu^{\mathrm{(ex)}}(n=0)=k_{\mathrm{B}}T\langle e^{\beta\varepsilon}|n=0\rangle$ is the contribution to the hydration free energy from solute interaction with the solvent outside the inner-hydration shell; $\varepsilon$ is the solute-solvent interaction energy, and the averaging $\langle\ldots|n=0\rangle$ is performed with solute-solvent thermally coupled, but with the inner-shell empty ($n=0$) of solvent molecules. As usual $\beta=1/k_{\mathrm{B}}T$, where $T$ is the temperature and $k_{\mathrm{B}}$ is Boltzmann’s constant.

Our earlier studies17, 18, 19, 20, 21, 22 establish $\lambda=5$ Å (Fig. 1) is a conservative definition of the inner shell. For $\lambda\approx 3$ Å the chemistry contribution is zero, uniquely identifying the domain excluded to the solvent. In our simulations, the inner-outer boundary is defined by a smooth, repulsive potential25, but this choice is inconsequential for the discussion below.

Previously, Hummer et al. 6 had carefully examined electrostatic contributions in the hydration of imidazole and imidazolium. We follow the same simulation setup ($NVT$ ensemble) and same potential model, but obtain the charging free energy using Eq. 1. The free energy of charging the imidazole is given by the difference of $\mu^{\mathrm{(ex)}}$ for imidazole and the analog with all charges set to zero ($\bf{Q=0}$). Since the packing contribution cancels in this difference, we need focus only on the chemistry term and the electrostatic contribution to the long-range term. Fig 2 (top panel) shows that the chemistry contribution is system size dependent, but the net electrostatic contribution is not. The excellent agreement with the value computed by Hummer et al. 6 confirms the internal consistency of our calculation.

Fig 2 (bottom panel) shows that the underlying chemistry (hydrophilic) and packing (hydrophobic) contributions are system size dependent, but these tend to compensate each other in the assessment of the net hydration free energy. However, the compensation is only partial. Thus, for example, for $N=216$, $\mu^{\mathrm{(ex)}}=-10.6\pm 0.2$, whereas for $N=512$, $\mu^{\mathrm{(ex)}}=-11.7\pm 0.2$. Note also that for $N=64$ packing dominates chemistry, a behavior found for hydrophobes26, 27: this artifact arises due to a constant volume constraint (supplemental information, SI). As our results show, this unphysical behavior is masked in the calculation of $\mu^{\mathrm{(ex)}}_{\mathrm{elec}}$.

The linearity of $\ln x_{0}$ and $\ln p_{0}$ with $v_{o}/L^{3}$ can be rationalized on the basis of two complementary perspectives. Analysis of fluctuations of number of solvent in a control volume (within a simulation cell) shows that a purely geometric factor of $(1-v_{o}/L^{3})$ factors from nonideal contributions 10, 11. This observation interpreted within a two-moment information theory model28, 29 then shows that $\ln x_{0}$ or $\ln p_{0}$ should have a leading order dependence on $v_{o}/L^{3}$. Alternatively, assuming ideal gas statistics gives the leading order dependence, we expect $\ln x_{0}$ (or $\ln p_{0}$) to depend as $\ln(1-v_{o}/L^{3})\approx-v_{0}/L^{3}$.

We next consider the hydration of the 56-residue protein $G_{B}$ (PDB: 2LHD), which is one member of a pair of conformational switch peptides. The protein is also net neutral, helping minimize corrections that are necessary in treating the hydration of a charged molecule 6. ($G_{B}$ is remarkable, for a single mutation flips the conformation of the protein from the $4\beta+\alpha$ fold to a $3\alpha$ fold 3.)

We first establish the reference $\mu^{\mathrm{(ex)}}_{\mathrm{elec}}$ (Fig. 3) for $G_{B}$ on the basis of standard coupling parameter integration (SI).

Within statistical uncertainties, $\mu^{\mathrm{(ex)}}_{\mathrm{elec}}$ including appropriate corrections (SI) is independent of the system size. We do not include the electrostatic finite size correction because the leading order monopole contribution, zero for $G_{B}$, only has a $R^{2}/L^{3}$ dependence 6, where $R$ is Born-radius of the solute. The results in Fig. 3 suggest that any finite size correction arising from the dipole contribution should be negligible. Further, just as we found for imidazole, $\mu^{\mathrm{(ex)}}_{\mathrm{elec}}$ obtained using the quasichemical approach is in excellent agreement with the result based on coupling parameter integration. Having established this consistency, we next consider the individual packing and chemistry contributions.

Figure 4 (top panel) shows that the (unfavorable) packing or primitive hydrophobic contribution becomes weaker with increasing system size as does the (favorable) chemistry contribution (Fig. 4, bottom panel). These trends are similar to that for imidazole.

Figure 5 gives the hydration free energy of the $G_{B}$ and its $\bf{Q=0}$ analog. The hydration free energy includes the long-range contribution, $\mu^{\mathrm{(ex)}}(n=0)$, which are only weakly dependent on the system size (SI).

From the extrapolated $L\rightarrow\infty$ value, we estimate the charging free energy $-628-33=-661\pm 3\,(2\sigma)$ kcal/mol noted in Fig. 3.

The results above definitively establish that the first hydration shell occupancy, and hence the hydration thermodynamics itself, is sensitive to the size of the simulation system. The system size dependence noted here is quite general and arises because the first shell is an open system that exchanges solvent with the bath. This feature, one that is made explicit by quasichemical theory 12, 13, 14, 15, has implications for all simulation studies of solvation, in general, and hydration, in particular.

The present work suggests that the observations by Karplus, Meuwly, and coworkers1 merits further study, but their rationalization of enhanced hydrophobic hydration in larger systems should be reconsidered. We find that the hydrophobic contribution is weakened with increasing system size, consistent with the intuition of enhanced solvent density fluctuations at the scale of the observation volume in a larger system. However, the hydrophilic (chemistry) contribution becomes less favorable as well. That is the solvent is less effective in prying apart or loosening the protein structure, and this effect dominates the weakening of hydrophobic hydration, as our recent studies emphasize 19, 20, 21, 22. This combined effect, which can be misinterpreted as a consequence of hydrophobic hydration22, should stabilize the conformation (assuming no change in intra-molecular interactions).

The criticisms by Gapsys and de Groot also merit serious consideration. They emphasize that the analysis in Ref. 1 cannot satisfactorily address the issue of system size effect, but they do not foreclose the possibility of system size effects in kinetics and thermodynamics of conformational transitions in that study. Importantly, for several model systems they find that the system size effect on the change in free energy along a conformational transition is small. We suspect that relative free energies of conformational transitions mask compensating effects of system size on hydrophobic and hydrophilic interactions noted here. Exploring this issue in biomolecular folding and assembly is left for future studies.

{acknowledgement}

We gratefully acknowledge computing support from National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract # DE-AC02-05CH11231. D. A. thanks Lawrence Pratt and Walter Chapman for helpful comments and encouragement. We thank Vyatautas Gapsys and Bert de Groot for helpful comments and clarifications.

{suppinfo}

Information on the simulation systems, quasichemical calculations, calculation of reference value of charging free energy, tabulated data, and analog of Fig. 2 for $G_{B}$.