All proteins contain groups capable of exchanging protons with their environment. We present here an approach, based on a rigorous thermodynamic cycle and the partition functions for energy levels characterizing protonation states of the associating proteins and their complex, to compute the electrostatic pH-dependent contribution to the free energy of protein-protein binding. The computed electrostatic binding free energies include the pH of the solution as the variable of state, mutual "polarization" of associating proteins reflected as changes in the distribution of their protonation states upon binding and fluctuations between available protonation states. The only fixed property of both proteins is the conformation; the structure of the monomers is kept in the same conformation as they have in the complex structure. As a reference, we use the electrostatic binding free energies obtained from the traditional Poisson-Boltzmann model, computed for a single macromolecular conformation fixed in a given protonation state, appropriate for given solution conditions. The new approach was tested for 12 protein-protein complexes. It is shown that explicit inclusion of protonation degrees of freedom might lead to a substantially different estimation of the electrostatic contribution to the binding free energy than that based on the traditional Poisson-Boltzmann model. This has important implications for the balancing of different contributions to the energetics of protein-protein binding and other related problems, for example, the choice of protein models for Brownian dynamics simulations of their association. Our procedure can be generalized to include conformational degrees of freedom by combining it with molecular dynamics simulations at constant pH. Unfortunately, in practice, a prohibitive factor is an enormous requirement for computer time and power. However, there may be some hope for solving this problem by combining existing constant pH molecular dynamics algorithms with so-called accelerated molecular dynamics algorithms.

It is well-known that the stability of proteins is pH-dependent because proteins contain groups capable of exchanging protons with their environment. In general, all conformational transitions of proteins might be expected to exhibit some pH-dependence. Using structural data for folded and partially unfolded forms of the Drosophila melanogaster engrailed homeodomain, we show that knowledge of the pK(a)'s of ionizable groups in different structural forms of a protein is not sufficient to determine pH-dependence of the free energy of its conformational transitions with accuracy better than 1 kcal/mol. We also describe a convenient computational approach for predicting expected errors when free energy of conformational transition is derived exclusively from the pK(a) values. This approach uses the full interaction free energy matrix between the ionizable sites, for the two forms of the protein, which can be determined when appropriate structural information is available.