Fitting multiple pKa's

I have three histidine residues that are next to each other in a protein. When I compute the fraction of unprotonated states for each individual histidine as a function of $$\mathrm{pH}$$ my titration curve has a non-Henderson–Hasselbalch (HH) shape. Hence I can't fit it using the HH equation for 1 site:

$$\frac{1}{1 + 10^{n(\mathrm{p}K_\mathrm{a} - \mathrm{pH})}}$$

in order to solve for the Hill coefficient $$n$$ and $$\mathrm{p}K_\mathrm{a}.$$

Obviously the three histidine $$\mathrm{p}K_\mathrm{a}$$'s are dependent on each other since they are very close to one another and are probably sharing the proton. So how do I go about fitting the curves? Should I sum the total unprotonated states in all three histidines and fit the curve with the three site equation

$$\frac{1}{1 + 10^{n_1(\mathrm{p}K_\mathrm{a1} - \mathrm{pH})}} + \frac{1}{1 + 10^{n_2(\mathrm{p}K_\mathrm{a2} - \mathrm{pH})}} + \frac{1}{1 + 10^{n_3(\mathrm{p}K_\mathrm{a3} - \mathrm{pH})}}?$$

When I use this equation I still don't get a good fit and I get unreasonable $$\mathrm{p}K_\mathrm{a}$$'s $$(-1.4, 7.7, -1.3).$$ Any suggestions?

As requested below, the matlab code I used for the fit is the following:

ft=fittype('(1/(1+10.^(a1*(p1-x))))+(1/(1+10.^(a2*(p2-x))))+(1/(1+10.^(a3*(p3-x))))', 'independent','x','dependent','y');
opts = fitoptions('Method','NonlinearLeastSquares');
[fitresult, gof] = fit( xData, yData, ft, opts );
h = plot(fitresult,xData,yData)


I posted values for p1, p2, p3 which should report $$\mathrm{p}K_\mathrm{a1},$$ $$\mathrm{p}K_\mathrm{a2}$$ and $$\mathrm{p}K_\mathrm{a3}.$$

• I suggest you post the code for the matlab function, at least. The fitting parameters you reportedly obtained do not match the shape of the curve in the figure. – Buck Thorn Dec 13 '19 at 6:50
• I second the comment above - please provide more detail of what and how you are doing the fitting? Are you fitting to the n values and allowing them to vary freely? It would make more sense to fit directly to the concentrations of 3H, 2H, H and no H with the constraint that they must all sum to 1. I would also worry about the effect of very high and very low pH on your protein. Is it stable across that range? The behavior should be similar to what one gets with a triprotic acid such as citric acid. For example, look at the chemistry libretext entry on "titration of a polyprotic weak acid". – Andrew Dec 13 '19 at 14:07
• As pointed out in one of the answers, your equation as written should reach a maximum of three at high pH not one. You need to put a correction factor to get from three to a scale of 0 to 1. More generally, fitting to that exponential form is extremely difficult with non-linear regression, even for a single HH equation. Replacing those denominators with the relevant expressions in terms of Ka’s and [H] will make the fit more robust. Even with that, fitting six independent parameters accurately with that small dataset will be hard. I would try setting all Hill coeffs to 1 on first try. – Andrew Dec 13 '19 at 20:42
• Do you have additional measurements than the 13 blue points represent? I'm not a biochemist, but to determine out three $pK_a$ (and not only one) with one measurement per unit of pH appears to me as coarse data, especially with the jumps for pH ~5 (in comparison to pH ~4 and ~6) and for pH ~7 in reference to the interpolation in red. As a suggestion: What about three runs each with measurements in increments of 0.1 units of pH, ploting both a) fraction of unprotonated in function of pH and b) the first derivative (change of fraction of unprotonated over change of pH value? – Buttonwood Dec 13 '19 at 21:34

Given your hypothesis that the three $$\mathrm{p}K_\mathrm{a}$$'s are dependent on one another and that experimentally you do not seem to be able to identify which histidines are protonated and which are not at any given time, it would seem that the simplest model to use is that of a triprotic acid such as citric acid. The triple HH equation model you used is only appropriate for three completely independent $$\mathrm{p}K_\mathrm{a}$$'s (such that the $$\mathrm{p}K_\mathrm{a}$$ at site 1 is unchanged by the protonation state of site 2).

In the dependent $$\mathrm{p}K_\mathrm{a}$$ model, we have four species (vs the eight in your model): $$\ce{H3P},$$ $$\ce{H2P},$$ $$\ce{HP}$$ and $$\ce{P},$$ where "$$\ce{P}$$" is the protein. There are three "apparent $$\mathrm{p}K_\mathrm{a}$$'s" none of which represents a true $$\mathrm{p}K_\mathrm{a}$$ of a specific site. They instead correlate to the combination of possible reactions that give a change in the overall protonation state. For example, the first $$\mathrm{p}K_\mathrm{a}$$ represents the change in state from $$\ce{H3P}$$ to $$\ce{H2P},$$ regardless of which of the three sites loses a proton.

We can then write four equations starting with mass balance:

\begin{align} [\text{total protein}] &= [\ce{H3P}] + [\ce{H2P}] + [\ce{HP}] + [\ce{P}]\tag{1}\\ K_\mathrm{a1} &= \frac{[\ce{H}][\ce{H2P}]}{[\ce{H3P}]}\tag{2}\\ K_\mathrm{a2} &= \frac{[\ce{H}][\ce{HP}]}{[\ce{H2P}]}\tag{3}\\ K_\mathrm{a3} &= \frac{[\ce{H}][\ce{P}]}{[\ce{HP}]}\tag{4} \end{align}

We can also define the fraction of sites unprotonated as

$$\phi = \frac{3[\ce{P}]+2[\ce{HP}]+[\ce{H2P}]}{3[\text{total protein}]}\tag{5}$$

It's relatively straightforward to use the four above equations to get expressions for each of the species concentrations in terms of $$[\ce{H}]$$ and $$[\text{total protein}]$$ only. Substituting those expressions into the bottom equation will give you an equation to which to fit your data.

If you post your data in tabular form, others can play around with different fitting approaches and give feedback.

• I understand that it might be an improvement but I am not sure this is entirely appropriate because the protons are not quite equivalent (the microscopic pKas are not identical at any given pH). – Buck Thorn Dec 13 '19 at 21:54
• They don’t need to be equivalent. In fact, if they are completely independent and have widely spaced pka’s, the apparent pka’s will converge on the microscopic pka’s – Andrew Dec 13 '19 at 22:03
• I took a liberty to brush up formatting a little and unify notations such as "[total protein]" and "total protein concentration", HH and H-H etc. for clarity. Also, I changed $[\ce{H3}]$ to $[\ce{H3P}]$ in eqn. (1). – andselisk Dec 14 '19 at 7:49

It makes no sense to solve for different pKa´s. They are constants specific to the acid groups of Histidine. Furthermore, the HH-equation describes the fraction of (in your case) dissociated acid in one acid/base equilibria. So if you sum the three steps, your Y-axis for one thing, has to go from 0 to 3. Just off the top of my head, it seems a better solution to multiply the HH-equations: At a given pH, a fraction of His has dissociated the first proton, and out of this fraction another fraction has given of one more proton etc.

The data look as if you have simply done a titration of histidine, with a first pKa at around 6, and another one at around 8.5 - 9. This is exactly what is expected from the tabulated values for histidine : pKa = 6.04 and 9.09. Why are you not satisfied about these measurements ? Your results are in perfect agreement with predictions form the tabulated values. It is no use trying to fit your points with the strange analytical curve that you propose.

• I believe the 9.09 pKa you are referring to is for the amino terminus. Within a polypeptide chain, this proton is not acidic, as the amine is now part of an amide bond. – Andrew Dec 13 '19 at 14:10