[R] non-linear regression and root finding

[Troels Ring]

Dear friends - I have a function for the charge in a fluid (water) 
buffered with HEPES and otherwise only containing Na and Cl so that [Na] 
- [Cl] = SID (strong ion difference) goes from -1 mM to 1 mM. With known 
SID and total HEPES concentration I can calculate accurately the pH if I 
know 3 pK values for HEPES by finding the single root with uniroot

Now, the problem is that there is some disagreement in the literature 
what is the correct value for the 3 pKs. I know the 3 pK values have the 
relationship pK1 < pK2 <- pK3 and for the most common formulation of 
HEPES I know the charge on fully protonated species is 2. Hence I can 
generate a huge number of pK values from uniform distribution by taking 
sort(runif(3,-1,9) and make sure there are no ties and then run all 
those triplets of pK values to find pH values and thereby find the 
lowest deviation vis a vis measured values. That works but requires many 
triplets and is not satisfying. Hence I wonder if I could somehow have 
non linear regression to find the 3 pK values. Below is HEPESFUNC which 
delivers charge in the fluid for known pKs, HEPTOT and SID. Is it 
possible to have root-finding in the formula with nls? (I know the 
precision asked for is extreme but it has worked well in really many 
applications).

All best wishes

Troels Ring, MD
Aalborg, Denmark


HEPESFUNC <-
   function(H,SID,HEPTOT,pK1,pK2,pK3) {
     XX <- (H^3/(10^-pK3*10^-pK2*10^-pK1)+H^2/(10^-pK3*10^-pK2)+H/(10^-pK3))
     IV <- HEPTOT/(1+XX)
     I <- IV*H^3/(10^-pK3*10^-pK2*10^-pK1)
     II <- IV*H^2/(10^-pK3*10^-pK2)
     III <- IV*H/10^-pK3
     H - kw/H + SID + I*2 + II - IV
   }


HEPTOT <- 0.050
SID <- c(-seq(10,1)*1e-4,0,seq(1,10)*1e-4)

pHobs  <- c(4.63,4.68,4.72,4.77,4.83,4.9,4.96,5.04,5.12,5.21,5.3,
             5.39,5.48,5.55,5.63,5.69,5.74,5.8,5.85,5.89,5.93)

pK1 <- -1; pK2 <- 3; pK3 <- 7.55 # literature values
pK3 <- 7.6; pK2 <- 2.96 # values eye-balled to be better

pH <- c()
for (i in 1:length(SID)) {
   pH[i] <- -log10(uniroot(HEPESFUNC,c(1e-20,1),tol=1e-20,maxiter=1e4,
                        HEPTOT=HEPTOT,SID = SID[i],
                        pK1=pK1,pK2=pK2,pK3=pK3)$root)
}

plot(SID,pH)
points(SID,pHobs,col="red")