Charged binding calculations
Authored by Gabriel Rocklin
Last updated July 2 2014
Overview
Many interesting processes involve ions changing chemical environment, such as the transfer of an ion from the gas phase to the solvated phase, the protonation of a protein side chain, or the binding of a charged ligand. Free energy simulations of these processes can be done in two ways: (1) the same simulated system can be used for both chemical environments, such as when a charged ligand is "pulled" from a protein binding site into the bulk solvent, or (2) different simulated systems can be used for different chemical environments, such as when a charged ligand is alchemically decoupled from its own "free in solution" box and alchemically inserted/coupled into the protein binding site in a different simulated system. In both cases, the computed free energies can depend on details of the calculation setup, such as the boundary conditions applied at the edge of the explicit solvent system (periodic, implicit solvent, or vacuum), the size of the periodic unit cell, and even completely arbitrary details of the solvent molecular topology.
These dependences on setup parameters are referred to as "finite size effects" because they are caused by the finite and extremely small size of the simulated system/periodic unit cell relative to the macroscopic size of the experimental system. Finite size effects are basically computational artifacts because they contribute to the simulated free energy change but not to a bulk experimental measurements. It is thus desirable to remove these free energy contributions before comparing simulated free energy results with experimental results.
The conventional way to do this is to introduce correction terms to convert the raw simulation results (which are specific to the particular choices of boundary condition, system size, molecular topology, etc) into corrected results based on ideal physical conditions, meaning that the solute is infinitely diluted in a nonperiodic bulk of solvent. Once simulation results have been corrected to match these ideal conditions, different simulation methods can be compared to each other without confounding artifacts, and simulation results can be compared with experimental results. (Of course, bulk experiments are not exactly performed under true ideal conditions either, but the non-idealities in experiments have little to do with finite size artifacts in simulations (periodic or otherwise) so it is still best to correct simulation results for these artifacts prior to comparisons with experimental results.)
This page provides short guidance on how to remove finite size artifacts for several common simulation approaches. For full details see the associated references. Note that the only types of artifacts considered here are those relating to the net charge of the investigated species (for instance, the charged ligand). For now, this page only covers simulations performed under periodic boundary conditions.
Periodic boundary conditions alter the electric field in the unit cell
The key difference between periodic boundary conditions and nonperiodic boundary conditions is how "zero" potential is defined. Under nonperiodic boundary conditions, zero is defined as the potential infinitely far from the system of interest.
[math]\displaystyle{ \Phi_{Nonper} (r=\infin) = 0 }[/math]
However, under periodic boundary conditions, zero potential is defined as the average potential over the periodic unit cell.
[math]\displaystyle{ \left \langle \Phi_{Per}(r) \right \rangle = 0 }[/math]
This requirement of average zero potential means that any perturbation to the average potential in the unit cell will alter the potential felt at every point in the unit cell even if no external electric field is produced. For example, a neutral sphere with a point negative charge at the center and a uniform shell of positive charge on the outside produces no electric field outside itself. However, there is a net negative potential inside the sphere. If this sphere is placed in a periodic unit cell, the overall negative potential inside the sphere produces a compensating uniform positive potential across the rest of the unit cell in order to maintain the zero average potential in the periodic unit cell.
Scenario 1: Alchemical annihilation of a charge in a periodic system
In an alchemical annihilation calculation, a charge is removed from a system over a series of simulated lambda intermediates. In most cases the charge is re-introduced in a different system, such as when a ligand is alchemically removed from bulk solvent and inserted into a protein binding site. In these calculations, the charge annihilation free energies depend on the electrostatic potential at the location of the charge, as calculated in the simulation model under periodic boundary boundary conditions. When a uniform potential is added to the box to maintain the required average potential of zero, the charge of interest feels this added artifactual potential. This is how the requirement that [math]\displaystyle{ \left \langle \Phi_{Per}(r) \right \rangle = 0 }[/math] alters binding free energies. To correct the simulated binding affinities, we need to understand the factors that alter [math]\displaystyle{ \left \langle \Phi_{Per}(r) \right \rangle }[/math].
Three main factors alter [math]\displaystyle{ \left \langle \Phi_{Per}(r) \right \rangle }[/math] and thus introduce finite size artifacts:
(1) Net charges in the system. All charges change the average potential in the unit cell. As a result, a "compensating background charge" is introduced which depends on the size of the periodic unit cell.
(2) The number density of water molecules in the system. Each water molecule behaves like the sphere in the above example - a negative charge on the inside and a positive charge on the outside. The magnitude of the compensating potential depends on the average potential inside the water molecule (which depends on which specific water topology is used), the number of water molecules in the box, and the size of the box.
(3) Complex charge distributions and excluded volume. If the simulation system contains a protein, the protein's charge distribution can change the average potential of the unit cell. The excluded volume of a large solute can also alter the average potential in the unit cell by changing the effective dielectric over some region of space from high (solvent like) to low (protein like).
A final type of finite size artifact does has nothing to do [math]\displaystyle{ \left \langle \Phi_{Per}(r) \right \rangle }[/math]. Unlike an infinite nonperiodic system, a periodic unit cell has only a limited amount of solvent. This means that an ion in a periodic system is "undersolvated" relative to the infinite dilution case. A correction for this effect is simple to apply.
A complete method for correcting these artifacts in protein-ligand binding calculations is given in:
Though the paper is long, the discussion section includes a summary "recipe" for applying the corrections to protein-ligand binding calculations. Just consider, one afternoon lost but new knowledge and better technique gained!
The major paper figures are summarized in this short presentation here.
Scenario 2: Protein-ligand binding in the same box
A common approach for calculating protein-ligand binding is to compute the potential of mean force (PMF) of "pulling" the ligand out of the protein binding side using artificial restraining potentials. The ligand is considered "unbound" once it is "suitably far" from the protein as to be considered non-interacting. While it is often believed this approach does not introduce finite size artifacts (because the simulation system can remain net neutral), this is not actually the case.
For a simple explanation of why, consider the solvation free energy of a spherical ion in water. When calculated under periodic boundary conditions, finite size artifacts perturb the uncorrected solvation free energy of the ion, and raw calculations need to be corrected to determine the solvation free energy of the ion at infinite dilution. It is well known and mathematically proven that these artifacts depend on the radius of the ion. For the full mathematical discussion, see:
This work demonstrates that spherical ion is marginally stabilized in a periodic box compared to in an infinite nonperiodic system, and that an ion of larger radius is stabilized to a greater degree. (This discussion does not consider effect 2 in the "Scenario 1" section.) Why, physically speaking, should a larger ion be stabilized by periodicity to a larger degree than a smaller ion? After all, if we only consider undersolvation (the last of the finite size effects described in "Scenario 1", above), both the large ion and small ion are undersolvated to the same degree, because undersolvation is defined as the absence of polarized solvent beyond the edge of the periodic unit cell. The reason the ion of larger radius is stabilized to a greater degree is because the excluded volume of the larger ion alters the average potential in the unit cell to a greater extent than in the case of the smaller ion. Each ion's excluded volume perturbs [math]\displaystyle{ \left \langle \Phi_{Per}(r) \right \rangle }[/math] and introduces a corresponding uniform offsetting potential to enforce [math]\displaystyle{ \left \langle \Phi_{Per}(r) \right \rangle = 0 }[/math]. This offsetting potential interacts favorably with the ion. The ion of larger radius causes a greater initial perturbation to [math]\displaystyle{ \left \langle \Phi_{Per}(r) \right \rangle }[/math], a greater offsetting potential, and a more favorable artifactual periodicity-induced stabilization.
Now, what does this mean for PMF calculations? As a charged ligand is pulled away from a protein, it effective shrinks in radius, as the excluded volume of the protein no longer occupies the region surrounding the ligand and instead occupies some distant region of space where the ligand's electric field is already greatly diminished. So, one of the components of the PMF calculation is the difference in solvation free energy between a "large ion" (the protein ligand complex) and a small ion (the free ligand distant from the protein), which is perturbed by periodicity as described above.
The above analysis should hold true regardless of how the ligand moves from the protein binding site into the solvent. In other words, a PMF approach, where the ligand is restrained at a series of increasing distances from the protein, should give the same binding free energy as an alchemical approach in which the ligand is alchemically decoupled from the protein and re-inserted into the same box at a large distance away from the protein. Both of these approaches will be affected by finite size artifacts to the same extent. The analysis also remains the same if, instead of inserting the ligand far from the protein as it is decoupled from the protein, a counterion is inserted instead. The nature of the re-inserted charge does not affect how the calculated free energy is altered by finite size effects.
One estimation of a finite size artifacts for one-box protein-ligand binding calculations has been published here:
The authors determine the periodicity artifact to be 1.09 kcal/mol in a 8 x 8 x 16 nm box (Table 3). Not trivial for those desiring precise comparisons with experiment!
An extended discussion of periodicity artifacts in PMF calculations is also contained in: PH Hünenberger and JA McCammon, Ewald artifacts in computer simulations of ionic solvation and ion–ion interaction: A continuum electrostatics study. J. Chem. Phys. 110, 1856 (1999)