Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2802180

Formats

Article sections

Authors

Related links

J Phys Chem B. Author manuscript; available in PMC 2010 November 26.

Published in final edited form as:

PMCID: PMC2802180

NIHMSID: NIHMS157141

The publisher's final edited version of this article is available at J Phys Chem B

We describe a method for determining diffusion-controlled rate constants for protein-protein association that explicitly includes the solution pH. The method combines the transient-complex theory for computing electrostatically enhanced association rates with an approach based on a rigorous thermodynamic cycle and partition functions for energy levels characterizing protonation states of associating proteins and their complexes. To test our method we determine the pH-dependent kinetics of association of the HyHEL-5 antibody with its antigen hen egg white lysozyme. It was shown experimentally that their association rate constant depends on pH, increasing linearly in the pH range 6–8 and saturating or even exhibiting a flat maximum in the pH range 8–10. The presented methodology leads to a qualitative agreement with the experimental data. Our approach allows one to study diffusion controlled protein-protein association under different pH conditions by taking into account the ensembles of protonation states rather than just the most probable protonation state of each protein.

All proteins contain ionizable groups that exchange protons with their environment. In solution, proteins fluctuate between available protonation states whose distribution is specific to a given pH. Protein ionizable groups are involved in intra-protein, protein-solvent, protein-protein and protein-ligand interactions influencing protein solubility, folding and catalytic activity.^{1} A proton transfer between ionizable groups and the surrounding solvent is a process with a pH-dependent equilibrium. To quantitatively describe such protonation equilibria in proteins, pK_{a}s of ionizable groups are often used. The pK_{a} of a specific ionizable group in a single protein or in a protein-protein complex can be significantly shifted relative to the pK_{a} of the corresponding isolated group in solution.^{2}^{,}^{3} Also, shifts in pK_{a}s between native and denatured states of proteins or between unbound and complexed forms are a source of the pH dependence of macromolecular stability or binding.^{4} As the exchange of protons modifies protein charge density, protein electrostatic properties are strongly related to the pH. Therefore, taking into account the protonation equilibria in proteins is essential for modeling their pH-dependent properties and functioning.

Association of proteins with each other and with other biomolecules is an important step in many biological processes. pH-dependent protonation equilibria determine both local charge distributions and protein net charges influencing protein-protein interactions and the kinetics of binding. Uptake or release of protons by protein ionizable groups can accompany their association and cause the association rate constant to be dependent on pH.

Effects of pH on protein association have not been widely studied. Standard theoretical techniques used to estimate rates of protein-protein association (such as Brownian dynamics) only partially take into account the pH effects;^{5}^{,}^{6} associating proteins are assigned single protonation states which are the most probable ones at a given pH either for isolated molecules or for their complexes. The facts that each protein exists as an ensemble of protonation states (with a number of different protonation states readily accessible by thermal agitations if their energies are sufficiently close to the energy of the most probable state) and that the protonation preferences of each of the associating proteins can be changed upon binding are neglected. Additionally, the pH at which the association occurs is not explicitly included during computations.

It was previously shown that the electrostatically enhanced rate of a diffusion limited association can be modeled with reasonable accuracy by^{7}^{–}^{9}

$${k}_{a}={k}_{ao}\xb7{e}^{-\frac{<{E}_{el}{>}^{*}}{{k}_{{B}^{T}}}}$$

(1)

where k_{a} and k_{ao} are, respectively, the association rates in the presence and absence of electrostatic interactions, < *E _{el}* >* is the average electrostatic interaction energy in a transient-complex ensemble that separates the bound and unbound states of proteins and is computed with the Poisson-Boltzmann (PB) model,

Even though protonation equilibria are crucial for many biological processes, the amount of available experimental data of protein-protein association kinetics under varying pH is limited. We chose to test our approach on a pair of proteins: the anti-hen egg lysozyme monoclonal antibody (HyHEL-5) and its antigen, hen egg lysozyme (HEL) (Figure 1). The large number of ionizable groups in both proteins and extensive titrations of mutually influencing residues makes the chosen system an interesting and challenging test case. The kinetics of association of these proteins under various conditions was previously studied by stopped-flow fluorescence polarization.^{12} It was found that the association rate of HyHEL-5/HEL complex varies with pH in the range 6.0–10.0. We determined the association rate constant in the HyHEL-5/HEL system using both the new method and the standard one. We show that including the protonation degrees of freedom and pH while evaluating the electrostatic interaction energies leads to significantly better agreement with experiment than the standard methodology.

The stereospecific binding of two proteins can be described as a consecutive two step process. First, proteins (A and B) approach each other via translational and rotational diffusion and reach a region in configurational space where their separations and orientations are close to those observed in the final complex. Such configurations of proteins constitute the transient-complex ensemble (C*). Second, the conformational rearrangement of proteins leads to the stereospecific final complex (C):^{7}

$$A+B\underset{{k}_{-a}}{\overset{{k}_{a}}{\rightleftharpoons}}{C}^{*}\stackrel{{k}_{b}}{\to}C$$

(2)

The overall association rate is:

$${k}_{on}=\frac{{k}_{a}{k}_{b}}{{k}_{-a}+{k}_{b}}$$

(3)

If the process is diffusion-controlled, k_{b}>>k_{−a} and k_{on}k_{a} approximation holds. In the transient-complex theory, the transient-complex ensemble represents configurations that are as close to the native complex as possible without the need to consider short range interactions and conformational rearrangements. If the two proteins in the transient-complex are severely constrained in relative translation and rotation, the diffusion controlled rate constant to reach the transient complex, k_{a} is given with eq 1.

The transient-complex ensemble is located at the outer boundary of the native-complex interaction energy well. The procedure to generate and identify the transient-complex ensemble follows that described previously.^{8}^{,}^{9} Proteins, in a configuration as in the native complex, are treated as rigid bodies with only six degrees of freedom: three for relative translation and three for relative rotation. The larger of the two proteins, HyHEL-5 is fixed in the laboratory frame. HEL translational and rotational coordinates are defined relative to the native-complex configuration built from the X-ray structure. Translation is represented in spherical coordinates by a relative displacement vector (ρ, θ, ϕ). Its norm, ρ, measures the separation between the surfaces of the two proteins. Rotation is represented by three Euler angles (α, β, γ) of the HEL-fixed frame, relative to the laboratory frame. The α angle describes the rotation around a HEL-fixed unit vector, β and γ are polar and azimuthal angles of the HEL-fixed vector in the laboratory frame. Coordinates are sampled from ρ [ρ* _{o}*, ρ

We generate random configurations of the two proteins by sampling the global coordinates (ρ,θ,ϕ,α,β,γ) from uniform distributions, so that all configurations are equally probable. Only the non-overlapping configurations are accepted and stored for further analysis. The overlaps are detected by testing all atoms of HyHEL-5 against all atoms of HEL. The distances between three classes of atoms are considered: hydrogens, polar (oxygens and nitrogens) and non-polar (carbons and other atoms). We assume that an overlap between atoms within the same or different classes occurs when their distance in the generated configuration is smaller than the minimum distance of such a pair in the X-ray native structure.

Non-overlapping configurations and contacts (native and non-native) are used to locate the boundary of the native-complex interaction energy well, as previously described.^{8}^{,}^{9} We choose a number of distinct atom pairs belonging to the interface between proteins and these native pairs form a list used to define contacts. Each atom is assigned a contact radius computed as half of the distance from its contact partner in the native structure of the complex. Native contacts occur if for a given non-overlapping configuration the distance between a native pair of atoms is smaller than the sum of its native distance plus 3.5Å. Non-native contacts are defined when any of the two atoms not forming a native pair are within the distance smaller than the sum of their contact radii plus 2.5Å. All steps of the procedure leading to the transient-complex ensemble were performed with a home-made software.

As stated above, the transient-complex ensemble separates the bound and unbound states of proteins. Its location can be identified by measuring the translational and rotational freedom of proteins in non-overlapping configurations. This can be done using σ_{α}(*N*), the standard deviation of the α angle among all non-overlapping configurations for a contact level N,^{8}^{,}^{9} or Ξ, the difference between σ_{α}(*N*) and its average at all contact levels lower than N. The maximum of Ξ coincides with the transient complex. The level of contacts characteristic for the transient-complex ensemble can be therefore identified by examining Ξ as a function of N:^{8}^{,}^{9}

$$\mathrm{\Xi}(N)=<{\sigma}_{\alpha}(N\prime ){>}_{N\prime <N}-{\sigma}_{\alpha}(N)$$

(4)

In order to compute E_{el} that includes protonation degrees of freedom, we use the previously described thermodynamic cycle.^{13} The electrostatic interaction energy for each configuration of proteins is:

$${E}_{el}={G}_{el}^{o,C}-\left({G}_{el}^{o,A}+{G}_{el}^{o,B}\right)+\Delta {G}_{el}^{nn}$$

(5)

where ${G}_{el}^{o,A},{G}_{el}^{o,B}$ and ${G}_{el}^{o,C}$ are electrostatic free energies of isolated proteins and of their configuration taken from the transient-complex ensemble. These terms explicitly include the solution pH as described further on. The last term, $\Delta {G}_{el}^{nn}$, computed for molecules with all ionizable groups of proteins in their neutral states, is pH-independent

$$\Delta {G}_{el}^{nn}=\frac{1}{2}{\displaystyle \sum _{i=1}^{{N}_{A}+{N}_{B}}{q}_{i}{\varphi}_{i}}-\left(\frac{1}{2}{\displaystyle \sum _{i=1}^{{N}_{A}}{q}_{i}{\varphi}_{i}}+\frac{1}{2}{\displaystyle \sum _{i=1}^{{N}_{B}}{q}_{i}{\varphi}_{i}}\right).$$

(6)

$\Delta {G}_{el}^{nn}$ is evaluated using the traditional PB approach where the solute is represented as a low-dielectric region with N_{A}, N_{B} and N_{A}+N_{B} partial atomic charges located at positions of atomic nuclei. The solute is surrounded by a high-dielectric continuum solvent containing dissolved electrolyte. The distribution of mobile charges (ions) in solvent is described according to the Boltzmann law. Dielectric polarizability is included implicitly by dielectric constants of the media. In eq 6, ϕ* _{i}* is the electrostatic potential at the location of the charge

The electrostatic free energy of a molecule with M ionizable groups can be computed as:

$${G}_{el}^{o}=-RT{\displaystyle \sum _{m}^{{2}^{M}}exp\phantom{\rule{thinmathspace}{0ex}}}\left(-\frac{{G}_{m}}{RT}\right)$$

(7)

where *G _{m}* is the free energy level associated with the m

$$\begin{array}{c}{G}_{m}\equiv G({x}_{1}^{m},\dots ,{x}_{M}^{m},pH,T)=2.303RT{\displaystyle \sum _{i=1}^{M}{x}_{i}^{m}\gamma i}\phantom{\rule{thinmathspace}{0ex}}\left(pH-p{K}_{a,i}^{mod}\right)\\ \hfill +{\displaystyle \sum _{i=1}^{M}{x}_{i}^{m}\Delta \Delta {G}_{ii}}+{\displaystyle \sum _{i=1}^{M-1}{\displaystyle \sum _{j=i+1}^{M}{x}_{i}^{m}{x}_{j}^{m}}\Delta \Delta {G}_{ij}^{int}}\end{array}$$

(8)

where ${\text{pK}}_{a,i}^{mod}$ is the negative decimal logarithm of the acidic dissociation constant of the isolated i^{th} ionizable group (known from experiments or quantum mechanical computations for small model compounds), ΔΔ*G _{ii}* is the free energy change for ionizing protein’s

It was shown before^{13} that it is possible to accurately compute the electrostatic free energy of a molecule with *M* ionizable groups without including all 2^{M} energy levels. Moreover, to determine the electrostatic free energy with accuracy better than $0.01\frac{kcal}{mol}$ only the energy levels from the lowest one (i. e., the most probable) G_{1} up to RT·ln(2)·M (L states) above has to be known:

$${G}_{el}^{o}={G}_{1}-RT{\displaystyle \sum _{m=1}^{L}\text{exp}}\left(-\frac{{G}_{m}-{G}_{1}}{RT}\right)$$

(9)

Computing of ${G}_{el}^{o}$ at a given pH, separately for HEL, HyHEL-5 and their transient-complex configurations, together with $\Delta {G}_{el}^{nn}$ allows one to determine pH-dependent electrostatic interaction energy defined by eq 5.

For comparison we also considered a "traditional" or standard approach where one calculates E_{el} for each configuration of the transient-complex ensemble assigning protonation states consistent with the lowest energy state determined for selected pH and for the X-ray structure of the complex:

$${E}_{el}^{trad}={G}_{el,1,X}^{o,C}-\left({G}_{el,1,X}^{o,A}+{G}_{el,1,X}^{o,B}\right)+\Delta {G}_{el}^{nn}$$

(10)

where indices 1 and *X* denote that only the lowest protonation state of the X-ray native complex is considered.

The coordinates of the HyHEL-5/HEL complex (Figure 1) determined at 1.7Å resolution were taken from the Protein Data Bank (accession code 1YQV^{17}). To obtain fully protonated molecules hydrogen atoms were added and their positions were optimized using the GROMACS package.^{18} Partial charges and radii were assigned to atoms according to the OPLS-AA force field.^{19}

Total number of 3.8·10^{7} random configurations of HyHEL-5 and HEL was generated. All these configurations were tested for clashes and 417894 non-overlapping ones were accumulated. By analyzing contacts in these non-overlapping configurations the transient-complex ensemble was characterized as described above.

Terms ΔΔ*G _{ii}* and $\Delta \Delta {G}_{ij}^{int}$ (eq 8) for HEL, HyHEL-5 and each of their transient-complex configurations were evaluated using MEAD (Macroscopic Electrostatics with Atomic Detail).

$\Delta {G}_{el}^{nn}$ and ${\mathrm{E}}_{el}^{trad}$ for each of the transient-complex configurations were computed with APBS^{22} by solving the PB equation on the 225Å×225Å×289Å grid with final spacing of 0.48Å. All electrostatic computations were performed at 150mM ionic strength and temperature of 298.15K. Dielectric constants were set to 4 for proteins and 78.54 for the solvent (water). Solute-solvent dielectric boundary was defined as the van der Waals surface. It was previously shown that for oppositely charged molecules the PB calculations with the van der Waals surface definition of a dielectric boundary give negative electrostatic interaction energies but with the solvent accessible surface definition the electrostatic interaction energies can be positive.^{23}

The transient complex is located at the outer boundary of the bound state. As described in the Methods section a unique definition of the transient-complex ensemble is possible using σ_{α} (the standard deviation of the rotation angle α) that measures the degree of rotational freedom in the set of generated non-overlapping configurations. Figure 2 plots the number of contacts (N) registered in 417894 configurations of HyHEL-5/HEL that passed the clash test versus σ_{α} computed for a given contact level. The Ξ function defined in eq 4, also presented in Figure 2, shows a maximum for 39 contacts.

The standard deviation of α, σ_{α} (continuous line), as a function of the number of contacts (N) determined for non-overlapping configurations of HyHEL-5/HEL. The bound and unbound states are characterized by a large and small number **...**

The number of configurations from the non-overlapping set with the contact level of 39 was 351; they were classified as the transient-complex ensemble. Representative examples of the transient-complex configurations are shown in Figure 3.

Representative configurations in the transient complex of HyHEL-5 (shown as van der Waals spheres) and HEL (cartoon models). Drawing was done using CHIMERA.^{11}

Scatter plots with the number of contacts as a function of relative translation ρ and rotation angle α prepared based on all random HyHEL-5/HEL configurations that passed the clash test are presented in Figure 4. For configurations characterized by 39 contacts, the relative translation ρ ranges from 3.7 to 6.6Å with an average of 4.6Å. The standard deviation of ρ observed in the transient-complex configurations is 0.5Å. The rotation angle α characterizing the transient-complex ensemble falls between −45° and 45°, however, we also observe weakly populated configurations with α ~ −180° and α ~ 180°. The N(α) plot is asymmetric mostly around α ~ −90° and α ~ 90°. Additionally, configurations with α around ±90° are characterized by smaller number of contacts. The average α computed over the transient-complex configurations is −1.2° (with a standard deviation of 9.6°) indicating that positive and negative α values are nearly equally sampled.

The total number of ionizable groups is 32 in HEL and 101 in HyHEL-5. Both our and previous calculations of residual pK_{a}s in HEL, HyHEL-5, and the HyHEL-5/HEL complex^{24} show that pK_{a}s of several residues located at the binding interface, including those involved in the Arg-Glu salt links, change significantly upon association. In agreement with^{24} we find that most acidic residues titrate in the pH range below 5.0 while most basic residues in the pH range above 9.0.

For pH between 5.0 and 10.0 we computed the net charges of molecules as an average over the ensemble of protonation states. The net charge of isolated HEL varies between +10e and +3e, that of isolated HyHEL-5 between +10e and −15e, and of HyHEL-5/HEL complex (as in the X-ray structure) between +18e and −13e (Figure 5). The resulting isoelectric points, pI, (pH of zero net charge) are ~11 for HEL (11.2 determined experimentally^{25}), ~7 for HyHEL-5 (consistent with pI of mouse immunoglobulines usually between 6.0 and 7.5^{26}^{,}^{27}) and ~9 for their complex.

Average net charges of HEL, HyHEL-5 and their complex (X-ray structure) computed as a function of solution pH.

At pH above 7.0, HEL and HyHEL-5 are oppositely charged. Experiments show that the maximal binding rate is observed for pH~8.0.^{24} At pH below 7 both molecules are positively charged but binding is still observed because local interactions between positively and negatively charged regions of both molecules and screening of interactions due to the non-zero ionic strength overcome the overall repulsion resulting from their positive net charges.

According to eq 9, the electrostatic free energy of a protein with multiple ionizable groups can be described as a sum of its free energy in the lowest protonation state and contributions from all other protonation states accessible at a given temperature. The latter contributions decrease the electrostatic free energy below its value in the lowest state G_{1} (see eq 9). If the number of easily accessible states is significant, the drop in the free energy below the lowest protonation state becomes substantial.

Figure 6 presents electrostatic free energy differences between the lowest protonation state G_{1} and ${G}_{el}^{o}$ of eq 9 computed for HEL and HyHEL-5 in the pH range 5–10. In general, the free energy differences observed for HEL are smaller than for HyHEL-5 what can be explained by smaller number of ionizable groups in HEL (32) than in HyHEL-5 (101). Large electrostatic free energy differences are found at low ($2.6\frac{kcal}{mol}$ at pH 5) and high pH ($3.9\frac{kcal}{mol}$ at pH 10) for the HyHEL-5 system and at high pH for HEL ($1.3\frac{kcal}{mol}$ at pH 10). Acidic ionizable groups titrate below pH 5 and basic groups start releasing protons at pH above 9; as a consequence the number of thermally accessible protonation states is raised at low and high pH. The number of accessible protonation states in the pH between 6 and 8 is lower because protonation/deprotonation of acidic/basic groups is less probable. This accounts for smaller decrease in the electrostatic free energy. For example, the number of states (*L* in eq 9) needed to evaluate the electrostatic free energy of HEL with accuracy better than $0.01\frac{kcal}{mol}$ is 258, 142 and 1157 for pH 5, 7, and 10, respectively. For HyHEL-5 the same accuracy is obtained with 68846, 21110 and 185897 states, respectively. Upon association the number of states accessible for both proteins and their distributions are likely to change. The number of states included in calculations of electrostatic free energies of transient-complex configurations varies (depending on a particular configuration and pH) from about 28000 (at moderate pH) to almost 500000 (at pH 10).

The difference between the energy of the lowest (i.e., the most probable) protonation state, G_{1}, and the electrostatic free energy, ${\mathrm{G}}_{el}^{o}$, computed for HEL and HyHEL-5 as a function of pH.

For transient-complex configurations, we present the distributions of electrostatic interaction energies, *E _{el}* (eq 5) and ${\mathrm{E}}_{el}^{trad}$ (eq 10) at different pH (see Figure 7). In general, the

Electrostatic interaction energy distributions computed for transient-complex configurations at various pH. Top plot: including the protonation degrees of freedom (E_{el}), bottom plot: after assigning a single protonation state (i.e., the most probable **...**

According to experiments^{12} the association rate constant for the interaction of HyHEL-5 with HEL varies with pH in the range 6.0–10.0 (see Figure 9). The experimentally measured k_{a} linearly increases with pH with a threefold rise in magnitude in the pH between 6.0 and 8.0. Above pH 8.0 k_{a} decreases but experiments were performed only at pH 10 so it is possible that the measured association rate exhibits a maximum between pH 8 and 10.

The association rate of HEL with HyHEL-5 computed as a function of pH. Experimental values are taken from;^{12} measurements were conducted at 150mM ionic strength and temperature of 298.15K).

To compare our calculations with experimental data, the basal rate for association, k_{ao} (eq 1), i.e., the association rate in the absence of electrostatic interactions, needs to be determined. Although it is possible to estimate the basal rate theoretically by simulating the transiational and rotational movement of proteins in the absence of electrostatic interactions (as described elsewhere^{7}^{–}^{9}), we treated k_{ao} as an adjustable parameter. We determined its value by finding a multiplicative factor that gives the best fit of the computed Boltzmann factors *e*^{−β<Eel>*}) and ${e}^{-\beta <{E}_{el}^{trad}{>}^{*}}$ (Figure 8) to the experimental curve in the pH range 6.0–8.0. The resulting basal rates are $2.74{\xb710}^{5}\pm {0.80\xb710}^{5}\frac{1}{M\xb7s}$ (based on fitting of ${e}^{-\beta <{E}_{el}^{trad}{>}^{*}}$) and $3.19\xb7{10}^{5}\pm 0.35\xb7{10}^{5}\frac{1}{M\xb7s}$ (based on fitting of *e*^{−β<Eel>*}). Both values fall in the ${10}^{5}-{10}^{6}\frac{1}{M\xb7s}$ range that is appropriate for stereospecific protein-protein association.^{28}^{,}^{29}

As shown in Figure 9, the traditional approach fails to reproduce the experimentally observed trend in the association kinetics. The slope of the computed curve in the pH region 6.0–8.0 is too small compared to the slope of the experimental curve. Moreover, the association rate constant computed with the traditional approach increases almost linearly in the whole studied pH range.

A much better result is obtained when protonation degrees of freedom and solution pH are explicitly considered. However, while in the pH range 6.0–9.0 the fit to experimental data is almost perfect, the drop in the association rate observed above pH 9.0 is significantly larger than experimental one. This may be due to a number of reasons. The presented approach neglects intra-protein conformational degrees of freedom. While protonation equilibria are treated explicitly and rigorously, the fact that proteins are frozen in their conformations found in the X-ray structure of the complex affects the results. Protein conformational freedom can be important particularly at low and high pH because titration of acidic or basic residues might be accompanied by local changes in the protein structure. Moreover, even some local changes in protein conformation (e.g. the side-chain reorientation) affect the computed energies of different protonation states.

Other limitations arise from the PB method. The dielectric model and the choice of solute dielectric constant may influence the computed shifts in the pK_{a}s of groups in the protein relative to their pK_{a}s in the corresponding model compound.^{24} Nevertheless, different dielectric models and dielectric constants are usually justified as a way to indirectly consider the internal flexibility of molecules.

We presented a computational approach to study kinetics of protein-protein association that explicitly takes into account proteins protonation equilibria and solution pH. We combined the transient-complex theory for protein association^{8}^{,}^{9}^{,}^{29} with explicit evaluation of electrostatic free energies of molecules with multiple protonation equilibria.^{13} The method was applied to determine the association kinetics of the HyHEL-5/HEL complex as a function of pH. Our results corroborate reasonably well with experimental data.^{12} The standard approach that we also tested fails to reproduce the pH-dependent trend in the association kinetics. While the standard approach gives almost monotonically increasing dependence of the association rate in the studied pH range, a maximum at pH~8 is observed in experiment.

The presented methodology is not free of limitations. Associating proteins are fixed in their crystal structure configurations while even local changes in protein conformation can affect the computed interaction energies. However, other methods widely used to study protein-protein association such as Brownian dynamics or PARE (Predicting Association Rate Enhancement)^{30} also suffer from this limitation. Moreover, these methods employ a less rigorous treatment of electrostatic interactions (as compared to the PB method), i.e., test charges or effective charge models and the Debye-Hückel approximation^{30}^{,}^{31} what may lead to additional uncertainties in predictions of association rate constants.

As proteins are represented with single structures, the coupling between protonation and conformational degrees of freedom is neglected. It is possible to include conformational flexibility during calculation of pH-dependent protein properties using either continuum electrostatic models (for example the MCCE method of Alexov and Gunner^{32}) or empirical approaches (the PROPKA method developed in the Jensen group^{33}). However, to incorporate conformational degrees of freedom in the presented methodology one would need to redefine the formalism used to compute the pH-dependent energy of interaction. Such modifications are underway.

We use the PB model to evaluate the electrostatic energies of protonation states. This model often overestimates the shifts in the pK_{a}s of ionizable groups in the protein relative to their pK_{a}s of the model compound. The results of the PB model may also depend on the force field used to assign atomic partial charges. While different methods exist to predict pK_{a}s and some of them outperform the PB approach,^{34} this is the only method that allows for rigorous evaluation of standard electrostatic free energies of molecules with multiple ionization sites.

The method used for calculation of free energy levels associated with different protonation states is based on the linear form of the PB equation. The electrostatic potential in the linear PB equation is a superposition of electrostatic potentials of each partial charge of the molecule. This property allows one to precalculate the interaction matrix,^{16} $\Delta \Delta {G}_{ij}^{int}$ (of the order equal to the number of ionizable groups in the molecule) that is later used to evaluate free energy levels according to eq 8. While the fact that the linear PB equation is used might prevent applications of our method to highly charged systems, the usage of non-linear PB equation requires performing PB computations on-the-fly with Monte Carlo generation of energy levels; such a task is computationally prohibitive or even impossible.

Notwithstanding its limitations, our approach can predict the association kinetics, its sensitivity to solution pH and help to analyze, understand and rationalize experimental observations. By comparing the presented method with the traditional one we find that the latter approach is reliable in the pH region where protonation/deprotonation events have low probability (neutral pH in case of HyHEL-5 and HEL). However, in pH regions, where protonation/deprotonation is more likely to occur (acidic and basic pH ranges in the studied case) our approach provides a more physically sound description and leads to a quantitatively better agreement with the experiment. The next step would be to check predictive capabilities of our approach.

M.D. and J.T. acknowledge support from Foundation for Polish Science, University of Warsaw (115/30/E-343/BST1345/ICM/2008 and G31-4), Polish Ministry of Science and Higher Education (N N301 245236), Fogarty International Center (National Institutes of Health Research Grant no R03 TW07318). J.A. is supported by University of Warsaw (BST). M.D. would like to thank Huang-Xiang Zhou and Sanbo Qin for helpful discussions.

1. Pace CN, Grimsley GR, J.M. S. J. Biol. Chem. 2009;284(20):13285–13285. [PubMed]

2. Alexov E. Proteins: Struct., Funct., Bioinf. 2004;56:572–584. [PubMed]

3. Mason A, Jensen J. Proteins. 2008;71:81–91. [PubMed]

4. Honig B, Nicholls A. Science. 1995;268:1144–1149. [PubMed]

5. Ermakova E. Biophys. Chem. 2007;130:26–31. [PubMed]

6. Ermakova E. J. Mol. Model. 2005;12:34–41. [PubMed]

7. Vijayakumar M, Wong K-Y, Schreiber G, Fersht A, Szabo A, Zhou H-X. J. Mol. Biol. 1998;278:1015–1024. [PubMed]

8. Alsallaq R, Zhou H-X. Proteins: Struct. Func. Bioinf. 2008;71:320–335. [PMC free article] [PubMed]

9. Qin S, Zhou H-X. J. Phys. Chem. B. 2008;112:5955–5960. [PMC free article] [PubMed]

10. Fogolari F, Brigo A, Molinari H. J. Mol. Recognit. 2002;15(6):377–392. [PubMed]

11. Pettersen E, Goddard T, Huang C, Couch G, Greenblatt D, Meng E, Ferrin T. J. Comput. Chem. 2004;25(13):1605–1612. [PubMed]

12. Xavier KA, Willson R. Biophys. J. 1998;74:2036–2045. [PubMed]

13. Piłat Z, Antosiewicz J. J. Phys. Chem. B. 2008;112:15074–15085. [PubMed]

14. Yang A-S, Gunner M, Sampogna R, Sharp K, Honig B. Proteins: Struct. Func. Genet. 1993;15(3):252–265. [PubMed]

15. Bashford D, Karplus M. Biochemistry. 1990;29:10219–10225. [PubMed]

16. Antosiewicz J, Briggs J, Elcock A, Gilson M, McCammon J. J. Comput. Chem. 1996;17:1633–1644.

17. Cohen G, Silverton E, Padlan E, Dyda F, Wibbenmeyer J, Willson R, Davies D. Acta Crystallogr., Sect. D. 2005;61:628–633. [PubMed]

18. Hess A, Kutzner C, van der Spoel D, Lindahl E. J. Chem. Theory Comput. 2008;4(3):435–447.

19. Kaminski G, Friesner R, Tirado-Rives J, Jorgensen W. J. Phys. Chem. B. 2001;105:6474–6487.

20. Antosiewicz J. Biophys. J. 1995;69(4):1344–1354. [PubMed]

21. Gilson M. Proteins: Struct. Func. Gen. 1993;15:266–282. [PubMed]

22. Baker N, Sept D, Joseph S, Holst M, McCammon J. Proc. Natl. Sci. USA. 2001;98:10037–10041. [PubMed]

23. Dong F, Zhou H-X. Proteins. 2006;65:87–102. [PubMed]

24. McDonald M, Willson R, McCammon J. Protein Eng. 1995;8(9):915–924. [PubMed]

25. Tanford C, Roxby R. Biochemistry. 1972;11:2192–2198. [PubMed]

26. Hamilton RG, M. R, Reimer C, Rodkey L. Electrophoresis. 1987;8:127–134.

27. Raman CS, Jemmerson R, Nall B, Allen M. Biochemistry. 1992;31:10370–10379. [PubMed]

28. Northrup S, Erickson H. Proc. Natl. Acad. Sci. USA. 1992;89:3338–3342. [PubMed]

29. Zhou H-X. Biophys. J. 1997;73:2441–2445. [PubMed]

30. Selzer T, Schreiber G. J. Mol. Biol. 1999;287:409–419. [PubMed]

31. Gabdoulline R, R.C. W. J. Phys. Chem. 1996;100:3868–3878.

32. Alexov G, Gunner M. Biophys. J. 1997;72:2075–2093. [PubMed]

33. Li H, Robertson A, J.H. J. Proteins. 2005;61:704–721. [PubMed]

34. Davies M, Toseland C, Moss D, Flower D. BMC Biochemistry. 2006;7:18–30. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |