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")
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
