Quantifying Protein-Protein Interactions in Molecular Simulations

Abstract

Interactions among proteins, nucleic acids, and other macromolecules are essential for their biological functions and shape the physicochemcial properties of the crowded environments inside living cells. Binding interactions are commonly quantified by dissociation constants K_d, and both binding and nonbinding interactions are quantified by second osmotic virial coefficients B₂. As a measure of nonspecific binding and stickiness, B₂ is receiving renewed attention in the context of so-called liquid-liquid phase separation in protein and nucleic acid solutions. We show that K_d is fully determined by B₂ and the fraction of the dimer observed in molecular simulations of two proteins in a box. We derive two methods to calculate B₂. From molecular dynamics or Monte Carlo simulations using implicit solvents, we can determine B₂ from insertion and removal energies by applying Bennett's acceptance ratio (BAR) method or the (binless) weighted histogram analysis method (WHAM). From simulations using implicit or explicit solvents, one can estimate B₂ from the probability that the two molecules are within a volume large enough to cover their range of interactions. We validate these methods for coarse-grained Monte Carlo simulations of three weakly binding proteins. Our estimates for K_d and B₂ allow us to separate out the contributions of nonbinding interactions to B₂. Comparison of calculated and measured values of K_d and B₂ can be used to (re-)parameterize and improve molecular force fields by calibrating specific affinities, overall stickiness, and nonbinding interactions. The accuracy and efficiency of K_d and B₂ calculations make them well suited for high-throughput studies of large interactomes.

Figure 1

Calculating the dissociation constant K_d and the second osmotic virial coefficient B_ij from simulations of two proteins in a box of volume V. The red protein has a single specific wedge-shaped binding site for the triangular blue protein. The light-blue protein configurations illustrate different interaction modes of the two proteins considered in the derivation of eqs 1 and 28. To obtain a K_d estimate independent of box size, we analytically extend the two-particle partition function for the simulation box by the contributions of an extension volume ΔV (gray shaded area) and perform the limit ΔV → ∞. We calculate B_ij from the probability p_v(V) that the two proteins are within a subvolume v (green), which is at least large enough to cover all protein–protein interactions (yellow shaded area).

Figure 2

Comparison of the accuracy of the insertion/removal method (ins/rem, black, solid lines) and the subvolume method (subvol, red, dashed lines) to estimate K_d and B_ij for three different protein pairs (top to bottom). The most likely estimates are indicated by horizontal and vertical dashed lines. The contour lines indicate the limits of the 25, 50, 75, and 95% confidence regions. The insertion/removal method (eqs 14 and 1 and the two-particle partition function from WHAM (Section 2.6.1), black) and the subvolume method (eqs 26 and 1, red) agree excellently with each other, and they have similar uncertainties. For UDG/Ugi, contour lines collapse on to a single line due to the strong correlation between the estimates for K_d and B_ij.

Figure 3

Box-size dependence of the binding probability p_b(V) is determined by B_ij and K_d via eq 15. We show simulation results (blue) for three protein pairs (top to bottom). Error bars indicate the blocked standard errors of the mean. The lines are predictions using eq 15 and estimates for K_d and B_ij obtained at a box volume = 3375 nm³ (magenta vertical line) using the insertion/removal method (black, solid lines) and the subvolume method (red, dashed lines).

Figure 4

Box-size dependence of the subvolume probability p_v(V) is determined by B_ij via eq 27. We show simulation results (blue) for three protein pairs (top to bottom). Error bars indicate the blocked standard errors of the mean. The lines are predictions using eq 27 and estimates of B_ij obtained at a box volume = 3375 nm³ (magenta vertical line) using the insertion/removal method (black, solid lines) and the subvolume method (red, dashed lines).

Figure 5

Contributions of binding and nonbinding interactions to B_ij = B_ij^(b) + B_ij^(u) for three protein pairs. We show estimates from the insertion/removal method in color and estimates from the subvolume method using larger symbols in gray. B_ij of the strongest binders is dominated by contributions of binding B_ij^(b) = −1/(2N_AK_d) such that the ratio of |B_ij^(b)/B_ij| is close to one (top). In these cases, nonbinding contributions to B_ij are relatively small, i.e., |B_ij^(u)/B_ij| ≪ 1 (center).

Figure 6

Finite-size correction gives box-size-independent dissociation constants K_d (eq 1, red symbols). The naive estimate of K_d (eq 8, blue symbols) suffers from finite-size effects and converges to the true value (gray horizontal line) for increasing box size. We illustrate this convergence by evaluating eq 8 for the predictions for p_b(V) from the insertion/removal method (black solid line) and the subvolume method (red dashed line). Approximately corrected estimates (eq 21, eq 13 of de Jong et al., green symbols) suffer from finite-size effects and converge to the true value for increasing box volume. Error bars have been obtained by resampling.

Figure 7

Relative difference between the approximate estimates K_d^′ (eq 21, eq 13 of de Jong et al.) and the box-size-independent estimates K_d (eq 1) for the dissociation constant as shown in Figure 6 (discs) as functions of the inverse box volume 1/V. This difference is proportional to 1/V and to the contribution of unbound states to the second osmotic virial coefficient, B_ij^(u) (eq 41, lines).