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")

? Mon, 6 Nov 2023 17:53:49 +0100 Troels Ring <tring at gvdnet.dk> ?????:> 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?Sure. Just reformulate the problem in terms of a function that takes a vector of predictors (your independent variable SID) and the desired parameters (pK1, pK2, pK3) as separate arguments and then returns predicted values of the dependent variable (to compare against pHobs): kw <- 1e-14 # I'm assuming pHm <- Vectorize(function(SID, pK1, pK2, pK3) -log10(uniroot(HEPESFUNC,c(1e-20,1),tol=1e-20,maxiter=1e4, HEPTOT=HEPTOT,SID = SID, pK1=pK1,pK2=pK2,pK3=pK3)$root)) (Yes, Vectorize() doesn't make the function any faster, but I don't see an easy way to rewrite this function to make it truly vectorised.) Unfortunately, nls() seems to take a step somewhere where crossprod() of the Jacobian of the model function cannot be inverted and fails with the message "Singular gradient". I wish that R could have a more reliable built-in nonlinear least squares solver. (I could also be holding it wrong.) Meanwhile, we have excellent CRAN packages nlsr and minpack.lm: minpack.lm::nlsLM( pHobs ~ pHm(SID, pK1, pK2, pK3), data.frame(pHobs = pHobs, SID = SID), start = c(pK1 = pK1, pK2 = pK2, pK3 = pK3), # the following is also needed to avoid MINPACK failing to fit lower = rep(-1, 3), upper = rep(9, 3) ) # Nonlinear regression model # model: pHobs ~ pHm(SID, pK1, pK2, pK3) # data: data.frame(pHobs = pHobs, SID = SID) # pK1 pK2 pK3 # -1.000 2.966 7.606 # residual sum-of-squares: 0.001195 # # Number of iterations to convergence: 15 # Achieved convergence tolerance: 1.49e-08 (Unfortunately, your code seemed to have a lot of non-breakable spaces, which confuse R's parser and make it harder to copy&paste it.) I think that your function can also be presented as a degree-5 polynomial in H, so it should also be possible to use polyroot() to obtain your solutions in a more exact manner. -- Best regards, Ivan

Your script is missing something (in particular kw). I presume you are trying to estimate the pK values. You may have more success with package nlsr than nls(). nlsr::nlxb() tries to get the Jacobian of the model specified by a formula and do so by applying symbolic or automatic differentiation. The multi expression function would probably not work. (That is one advantage of nls(), but it has an unstabilized solver and unless carefully called will use a simple forward derivative approximation.) There is also nlsr::nlfb() that can use functions for the residuals and jacobian. Your function is messy, but likely the jacobian can be developed with a bit of work. Some of the symbolic derivative features of R can help here. Alternatively, there are several numerical approximations in nlsr and the vignette "Intro to nlsr" explains how to use these. Note that simple forward and backward approximations are generally not worth using, but central approximation is quite good. The solver in the nlsr package allows of bounds on the parameters. Since you want them ordered, you may want to transform pK1 < pK2 <- pK3 to pk1, deltapk2, deltapk3 so pK2 = pk1+deltapk2 and pk3 = pk2 + deltapk3 = pk1 + deltapk2 + deltapk3 and all values bounded below by 0 or possibly some small numbers to keep the paramters apart. I won't pretend any of this is trivial. It is VERY easy to make small errors that still allow for output that is disastrously incorrect. If it is possible to get a single "formula" as one expression even if spread over multiple lines, then nlxb() might be able to handle it. J Nash (maintainer of nlsr and optimx) On 2023-11-06 11:53, Troels Ring wrote: ....> 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")