On the vibrational free energy of hydrated proteins

In order to sample the configurational space, MD simulations were performed with Gromacs [25] version 4.6.3, using the Amber 99SB-ILDN forcefield [26] and the TIP3P water model [27]. Short range electrostatic and van der Waals interactions were cut off at 12 Å, long-range electrostatics being handled through the particle mesh Ewald method [28].

Each considered protein was embedded in a water box, its boundaries being at least 10 Å away from the protein. Sodium and chloride ions were added so as to neutralize the charge of the system and to reach a salt concentration of 150 mM L⁻¹, as often done [29–31].

The system was then relaxed using steepest descent minimization, with harmonic restraints (force constant of 1000 kJ mol⁻¹ nm⁻¹) on protein heavy atoms, until a maximum force threshold of 1000 kJ mol⁻¹ nm⁻¹ was reached. The solvent was next equilibrated at 300 K, first during 500 ps with a control of both volume and temperature, using Berendsen thermostat [32] and a coupling constant of 0.1 ps, then during 500 ps with a control of both pressure and temperature, using Parrinello–Rahman barostat [33] and a coupling constant of 2 ps. Finally, restraints were removed and 100 ns of simulation at 300 K were performed, with a timestep of 2 fs, all bonds being constrained with the LINCS algorithm [34].

From each simulation, 100 snapshots were picked (one per nanosecond), the n_clos closest water molecules were retained (n_clos ranging between 0 and 1000) and the corresponding hydrated conformers were energy-minimized with Tinker [35] version 6.2, using the Amber 99SB forcefield, L-BFGS minimization being performed until a gradient of 0.001 kcal mole⁻¹ Å⁻¹ was reached (with the minimize program), followed by truncated Newton minimization down to a gradient of 0.000 01 kcal mole⁻¹ Å⁻¹ (with the newton program), such a strict convergence criterion being required in order to obtain the expected six zero-frequencies [36–38].

For each energy-minimized conformer, normal mode frequencies, ν_i being the ith one, were then obtained by diagonalizing the mass-weighted Hessian matrix (with the vibrate program), allowing to obtain A_vib, the vibrational free energy, which is such that [18]:

$$\begin{equation}{A}_{\text{vib}}=\sum ^{{n}_{\mathrm{dof}}}\frac{1}{2}h{\nu }_{i}+\sum ^{{n}_{\mathrm{dof}}}kT\enspace \mathrm{ln}\left(1-{e}^{-\frac{h{\nu }_{i}}{kT}}\right),\end{equation} \tag{ 1 }$$

where n_dof is the number of degrees of freedom (n_dof = 3N − 6, N being the number of atoms of the system), h and k, respectively, the Planck and Boltzmann constants, T, the temperature, being set to 300 K.

In several previous studies (e.g. [19, 21, 22, 39]), only the contribution of the vibrational entropy to A_vib was considered. It is given by [18]:

$$\begin{equation}-T{S}_{\text{vib}}=\sum ^{{n}_{\mathrm{dof}}}kT\enspace \mathrm{ln}\left(1-{e}^{-\frac{h{\nu }_{i}}{kT}}\right)-\sum ^{{n}_{\text{dof}}}\frac{h{\nu }_{i}}{{e}^{\frac{h{\nu }_{i}}{kT}}-1}\end{equation} \tag{ 2 }$$

while the vibrational enthalpy is a follows:

$$\begin{equation}{H}_{\text{vib}}=\sum ^{{n}_{\mathrm{dof}}}\frac{1}{2}h{\nu }_{i}+\sum ^{{n}_{\text{ddl}}}\frac{h{\nu }_{i}}{{e}^{\frac{h{\nu }_{i}}{kT}}-1}\end{equation} \tag{ 3 }$$

the first sum corresponding to the zero-point energy.

As shown in figure 1(a), when 100 conformers of the BPTI picked from a 100 ns MD trajectory (see section 2.1) are energy-minimized, the standard deviation of their vibrational free energy (equation (1)) is much lower (4.2 kcal mole⁻¹) than the standard deviation of their potential energy (18.9 kcal mole⁻¹). However, when the energy of each BPTI conformer is minimized together with the 1000 water molecules that are the closest to this conformer, that is, all those that are at most ≈8 Å away from a protein heavy atom, a significant anti-correlation (−0.79) between the vibrational free energy and the potential energy of the energy-minimized conformers is observed (see figure 1(b)), meaning that low-energy hydrated conformers of BPTI are on average significantly more rigid than high energy ones.

Zoom In Zoom Out Reset image size

In agreement with previous studies [18, 19, 22], and as illustrated in figure 2 by comparing the highest and lowest vibrational free energy conformers of BPTI, most of the variation of the vibrational free energy is due to modes with frequencies below 1000 cm⁻¹. Note however that in this case, while the contribution of the vibrational entropy (equation (2)) is indeed dominant, the contribution of the vibrational enthalpy (equation (3)) is far from being negligible (26% of the difference, in the case of the conformers hydrated with 1000 water molecules), coming mostly from modes with frequencies ranging between 500 cm⁻¹ and 1000 cm⁻¹.

Zoom In Zoom Out Reset image size

On the other hand, the zero-point energy (equation (3)) represents 48% of the vibrational free energy difference between the pair of hydrated conformers, in line with the hypothesis that this term is a key component of the energy of macromolecular systems [46, 47], like in the case of small molecules [48–50].

In order to assess the importance of the anti-correlation between the vibrational free energy and the potential energy of energy-minimized conformers, a least-square fit of the data was performed, by assuming that:

$$\begin{equation*}{A}_{\text{vib}}\approx a{V}_{\mathrm{min}}+b,\end{equation*}$$

where V_min is the minimized potential energy of the considered conformer, a being the slope of the vibrational free energy—potential energy relationship.

Together with the BPTI, four other proteins of various folds and sizes were studied (see figure 3), retaining the 0, 10, 50, 100, 150, 200, 300, 400, 500, 600, 700, 800, 900 or 1000 closest water molecules of each MD conformer for the analysis (see section 2.2).

Zoom In Zoom Out Reset image size

As shown in figure 4, like for BPTI (figure 1(a)), when no water molecule is retained, the slope of the vibrational free energy—potential energy relationship has a small value, ranging between -0.03, for BPTI, and -0.15, for the lysozyme. On the other hand, when there is more than 3.8 water molecules per amino-acid residue in the considered system, the slope is always below −0.3, confirming that low-energy conformers are on average more rigid than high-energy ones. Moreover, the fact that the value of the slope does not seem to vary significantly when there is more than 3.8 water molecules per residue suggests that this rigidity is mostly due to the interactions between the protein and the water molecules belonging to its hydration shell. However, as illustrated in figure 5 in the case of the lowest (7416 kcal mole⁻¹) and highest (7466 kcal mole⁻¹) vibrational free energy conformers of BPTI, when ≈3.8 water molecules per residue are considered, the network of hydrogen bonds established between the water molecules all around the protein is almost continuous.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Interestingly, 3.8 water molecules per residue corresponds to 0.6 gram of water per gram of protein, that is, the order of magnitude of the amount of tightly bound water molecules found in the hydration shell of proteins [51] by various techniques like calorimetry [52] or microwave dielectric measurements [53, 54], in line with the idea that such water molecules are structural ones, that is, that they actually belong to the protein structure [55], even though their mean residence time can be quite short [56–58], namely, on the sub-nanosecond time-scale, in spite of a few notable exceptions [59–61].

In order to assess if the vibrational free energy could prove helpful for pinpointing physically relevant protein conformers, eight models proposed for target 624 during the 9th CASP experiment [62] were analyzed as above, together with the actual crystallographic structure (see figure 6). Target 624 was chosen because it proved possible to refine the two best models proposed (ts172-1 and ts172-4) through MD simulations on the sub-millisecond time-scale, the C_α root-mean-square (RMSD) to the native structure decreasing from ≈5 Å down to about 1 Å [8].

Zoom In Zoom Out Reset image size

The five models proposed by the group of David Baker (group 172) were included in the present analysis, as well as three models with an RMSD to the native structure of around 7 Å, proposed by three other groups (ts236-2, ts345-1, ts484-4). As shown in figure 7, except ts172-5, all proposed models remain relatively stable during a 100 ns MD simulation, with an RMSD to the initial structure ranging between 2 and 4 Å. Note however that the crystal structure is, according to this criterion, obviously the more stable one, with an RMSD staying below 1.5 Å during the whole simulation.

Zoom In Zoom Out Reset image size

As shown in figure 8, when no water molecule is retained (left column), several models have an average energy similar to, or even better than the average energy obtained starting from the crystal structure, namely: ⟨V_min⟩ = −0.01 kcal mole⁻¹, for ts172-1, ⟨V_min + A_vib⟩ = −7.5 kcal mole⁻¹, also for ts172-1, +3.3 kcal mole⁻¹, for ts345-1 or +6.1 kcal mole⁻¹, for ts172-2, the averages being obtained with 100 energy-minimized MD conformers (section 2.2), the zero corresponding to the average value obtained with the crystal structure.

Zoom In Zoom Out Reset image size

When water molecules close to the conformers are included in the analysis, the gap between the average energy of the proposed models and the average energy obtained starting from the crystal structure increases dramatically, being around 25 kcal mole⁻¹ when the vibrational free energy is considered on top of the potential energy, when n_clos ⩾ 200 (see figure 8(b)). However, when only the potential energy is taken into account, the gap is also large (see figure 8(a)), the smallest values being ⟨V_min⟩ = +9.5 kcal mole⁻¹, for ts484-4 (with n_clos = 400) and +14.1 kcal mole⁻¹, also for ts484-4 (with n_clos = 300).