Free-energy calculations

The standard binding free energy of an ion to a protein is expressed as

where the first term is the free-energy change in translocating the ion from bulk water to the binding site and the second term represents the loss of translational entropy during this process. The last term is used when the protein residues coordinating the ligand are restrained to improve convergence in the calculation of ΔG_int. If no restraints are applied, ΔΔG_c,P = 0. Assuming a Gaussian distribution for the ion positions in the binding site, the translational free-energy difference is estimated using (33)

where V₀ = 1661 Å³ is the reference volume for the standard concentration of 1 M, e is Euler’s number, and the σ-values are the principal root mean-square fluctuations of the ion positions. The restraint free energy is calculated using the thermodynamic integration (TI) approach of Cecchini et al. (34):

Here, |X − X₀|² is the mean-squared displacement of the restrained residues from a reference structure at a given spring constant k. We find that a value of 1.0 kcal/mol/Ų per atom is sufficient to reduce the fluctuations of the coordinating residues. The TI calculations are performed using 11 windows, and each window is simulated 200 ps for equilibration and 800 ps for production. To obtain ΔΔG_c,P, two separate calculations are performed, one with the ion at the binding site and another with the ion in bulk, i.e., ΔΔG_c,P = ΔGc,Pion−bulk−ΔGc,Pion−site.

The free energy of translocation is calculated using both free-energy perturbation (FEP) and TI (35). The FEP calculation is performed first, and the trajectory is used as the starting point for the TI calculations. As in the restraint free energy, the free-energy change between the states with the ion at the site and in bulk are calculated, i.e., ΔΔGint=ΔGintion−bulk−ΔGintion−site. However, only one calculation is needed if we create and destroy the ion simultaneously in the same system. For the FEP calculations, we use 66 exponentially spaced λ values because this was shown in previous studies to be adequate for charged molecules (14). At each λ value, we perform 50 ps of equilibration, followed by a further 50 ps of production run. In the TI method, the ensemble of the free-energy derivative is calculated at each λ value. The free energy is recovered by integrating the free-energy derivative with respect to the λ values:

For charged molecules, the integral can be done using Gaussian quadrature. A seven-point quadrature is used because this was shown to be adequate in a previous study (36). The starting point for each TI window is obtained from the corresponding FEP calculation, and each window is run for a total of 1 ns, with the first 200 ps discarded as the equilibration phase. The convergence of the TI calculations for the ions is shown in Fig. S1. The errors in free energies are estimated from block data analysis.

For the substrate, the standard binding free energy is expressed as

which includes the free-energy loss due to rotational entropy ΔΔG_rot, and the conformational restraints are applied on the substrate instead of the protein, as indicated in the last term, ΔΔG_c,L. The rotational free-energy difference is calculated in a similar manner to Eq. 2 (33):

Here, σ_ϕ-values are the rotational root mean-square fluctuations of the substrate at the binding site calculated using the quaternion representation. For the conformational restraint, a spring constant of 5.0 kcal/mol/Ų is applied to the heavy atoms of the substrate. The bulk calculation is performed separately, with the substrate placed in a 25 Å cubic water box with a neutralizing ion. As in previous calculations (14, 18), the translocation free energy is split into three stages to improve convergence:

The first term on the right-hand side is the electrostatic contribution, and the last two terms are the LJ contributions. The LJ term is split into the backbone (LJ-bb) and the side-chain (LJ-sc) parts to prevent a water molecule getting trapped behind the substrate in the binding site, as noted previously (14). For the electrostatic free energies, we report the values obtained with the TI method with 0.5/1.0 ns for equilibration/production. For the LJ free energies, we report the FEP values instead because of the lack of convergence with quadrature windows. A total of 52 and 37 λ-values are used for the LJ-bb and LJ-sc calculations, respectively. For each window, we run 100/100 ps for equilibration/production (see Fig. S2 for the convergence of the forward and backward transformation for both electrostatic and LJ free energies). A soft-core potential is used in the LJ calculations with a shifting coefficient of 6.0 (37).

Umbrella sampling simulations are performed for the Na1′ → Na3 transition of Na⁺ to find the potential of mean force (PMF) and study the conformational changes caused by this transition. The reaction coordinate is chosen along the vector from the Na1′ site to the Na3 site, and umbrella windows are generated at 0.5 Å intervals using k = 10 kcal/mol/Ų for the harmonic biasing potential. Overlaps between the samples of neighboring windows are found to be >10%, which assures a robust construction of the PMF using the weighted histogram analysis method. Each window is simulated for 8 ns. Convergence of the PMF up to the barrier is obtained after 3 ns, so the final PMF is constructed from the last 5 ns.