Talanra, Vol. 32, No. 6, pp. 461-465,
1985
0039-9140/85 $3.00 + 0.00 Pergamon Press Ltd
Pnnted in Great Britain
MATHEMATICAL OPTIMIZATION IN THE DETERMINATION OF THE DISSOCIATION CONSTANTS OF A TRIBASIC ORGANIC ACID BY 13C NMR SPECTROSCOPY ANNE-MARIE Department
of Chemistry
and Biochemistry,
(Receitjed 6 September
SALONEN
University
1984. Revised I5 December
of Turku,
SF-20500
Turku,
Finland
1984. Accepted 11 January 1985)
Summary-The calculation of dissociation constants from the chemical shifts of ‘% NMR spectra leads to a complicated non-linear equation. Two different mathematical methods for solution of this equation have been chosen-an iterative step method and a matnx pseudo-inversion method. When the iteration method is used the initial guesses for the parameters, the initial value of the step size and the escalation of the iteration must be optimized. For comparison the matrix pseudo-inversion method was used because it gives a unique result. With the optimized step method the results were as accurate or even better than those obtained with the matrix method. Although it takes time to optimize the system, the step method is the more suitable method of solving the problem. The matrix inversion can be done only with a computer with I3 significant digits and exponent capacity larger than +38.
In the case of strongly
overlapping ionization proApKa < 2.7,’ it has been difficult with classical methods to determine the dissociation constants of polyprotic acids.2 Computers have proved their worth for resolving this problem. The purpose of this paper is to show that complicated non-linear least-squares minimization problems can be solved even with micro-computers. An organic tribasic acid, 6-(4-carboxy-2-oxabutyl)d-ethyl-4,8_dioxaundecanedioic acid, C,5H2609. with both ApK,, values -0.6, has been examined as an example.
cesses,i.e.,
FORMULATION
p =
I
V-k-,A’-IW+l’.i=l,...,n WA
&=K,K2...K,= fi K, ,=I The total concentration
Ca=[H,A]+[H,_,A-]+...+[A”-] = and from the expressions stants, [H,_,A’-]
&bs= Cd,
H,A$H+
+H,m,A-
H,A+2H+
+H,_,A*-
H,A tijH+
+ H,m,A’-
for the dissociation
(2) con-
. dissociation
‘=’ n 1+ c
,_ I
B,lW+l’
Use of these expressions easily gives a general expression for the mole fractions: [H+]n-’ fi K, ,=o
+A”-
cumulative
= $$
i: [K,A’-1 [H,A]=
‘=;O(IH+y-,fiK,y
The corresponding are
,+-,A’-1
By substitution of these concentrations in equation (2) and rearrangement, the concentration of undissociated acid is found to be
(1)
where for the ith species x, is the mole fraction and 6, the chemical shift relative to the standard.4 To calculate the mole fractions it is first necessary to derive the concentrations of the species. For an acid with n dissociable hydrogen ions the following equilibria exist:
H,Atini;+
of the acid is
OF THE PROBLEM
The chemical shifts in the NMR spectra of acids and bases are a function of PH.~ For systems involving several ionized species the observed chemical shift is
.
’
constants
J =”
“‘.”
where & is defined as unity for convenience in notation, and x0 is the mole fraction of undissociated acid. Substitution of the mole-fraction expressions into equation (1) gives
ANNE-MARIE SALONEN
462
and
Matrh pseudo-inversion method After rearrangement, the form
equation
(3) can be written
in
(4)
MATHEMATICAL TREATMENT EQUATION (4)
OF
Nb=h
Optimized iterative step method To solve equation (4) two programs have been written in FORTRAN, program NDISS for the cumulative dissociation constants and program KDISS for the corresponding stepwise constants. Subprogram ST60@ was chosen for iteration and subprogram QUEN06 for error calculation. The optimization has been done by changing the initial guesses of the parameters, the initial step size and the escalation* of the iteration. The optimization criterion is minimization of the sum of the squares of the differences between the observed and calculated chemical shifts:
F = ~&bs -
For comparison with the iteration method the matrix method was used because it gives a unique result. The program was written in BASIC and the computer used was a WANG 2200 (Wang Laboratories, Massachusetts) with 13 significant digits and exponent capacity f99. Equation (4) must be changed into the matrix form
where h is a column vector of ones, N is a k x n matrix, k is the number of parameters, n the number of data-points and b is a column vector of the parameters. Because N is not a square matrix, when n > k, N must be multiplied by its transpose, NT. The matrix NTN is then a square matrix and the solved form of equation (5) is b = (N’N) ’ NTb (6) where (NTN)-‘NT
The calculations were run with an APPLE II+ (Apple Computer Inc., California) equipped with a CP/M-FORTRAN translator in addition to the standard Applesoft BASIC interpreter. The FORTRAN translator uses 6 or 7 significant digits and the exponent capacity is less than k38.
J
(5) it is possible
to write
cov(b) = cov[(NTN)- ‘NTh] According
to the theorem
cov(AX) = Acov(X)A’ we can write cov(b) = (NTN))‘NTcov(b)[(NTN))‘NT]T
The standard deviations of the parameters estimated by Quenouille’s method,6 n-l T,F;
of N.
Error estimation9
Error estimation
s(B) -
is the pseudo-inverse’
For the matrix equation the covariance matrix
&IIC)*
(5)
n (a --J*;
a =k,$
= (NTN)-‘N’cov(b)N(N”N)‘.
were
5,
It is assumed that the error is normally distributed, with mean zero and variance 6’. and h, elements are independent with equal errors, so cov(h) =,fI
where n denotes the number of data-points, and zi, is the value of a parameter when function F has been minimized with omission of the ith data-point. To find the errors of the p/3- and pkl-values the variance’ of an arbitrary function, Z = $(a, b), must be used:
wherej’= sum of squares/degrees the identity matrix. Thus
of freedom
and I is
cov(b) = (NTN)NTjIN(NTN)-’ and
var(Z) = si
,f = @ -
WT(h - Nb) n-k
to yield var(pp)
=x$
=
log e ’ - ___ si B
( >
*Escalation refers to the multiplication factor used in the Iteration.
It is sometimes
called
the magnification.
The errors of the p/?- and pK-values have been calculated in the same way as in the iteration method. EXPERIMENTAL
The ‘jC NMR-titrations FX-60
FT NMR
were done at 298 k spectrometer operating at
1K on a Jeol 15.03 MHz,
Dissociatton
constants
of a tribasic
organic
acid
463
CH 2-0-CH,CH,COOH ‘CH,-0
Fig. 1. The structural
formula
-
CH,CH,COOH
of the acid C,sH,,O,.
using 8000 data-points with a precision of k0.06ppm. Aqueous samples were placed in a lo-mm probe inside which was a 5-mm probe of DzO + DSS as an external standard. Potenliometric system
PH
Values of pH were measured at 298 k 0.5 K with a Radiometer ABL-3 potentiometer equipped with an Orion model 91-03-00 combined semimicro glass electrode assembly, standardized with NBS buffer solutions (pH 4.008 and 6.865 at 298 K). Sample preparation The concentration of the acid was 0.0501M (k4 x 10-4M). The samples were delivered by micropipettes to give a total volume of 1.5 ml and the titrant used was a standard 0.9989M (+ 6 x lo-“M) CO,-free potassium hydroxide solution (Merck “Titrisol”). The volume change has no effect on the chemical shift.“’
RESULTS AND DISCUSSION
The structural formula of the acid is shown in Fig. 1 and the ‘)C NMR spectrum of the acid is shown in Fig. 2; numbering of the peaks refers to the scheme in Fig. 1 and the external standard peaks are marked with an x. The chemical shifts of the carboxyl groups-peak number 7-were the object of this study. The chemical shifts of species H,A and A3- were measured in highly buffered solutions and 6(H,A-) and 6(HA*-) were calculated by the second method devised by Loewenstein and Roberts.” The 13CNMR titration curve, 6 us. pH, is plotted in Fig. 3. In Tables 1-3 it is quite easy to find the optimal values for the calculations of the dissociation constants: they are
Fig. 3. ‘jC NMR
titration
7
spectrum
chemical
shift US. pH.
denoted by asterisks. The optimum values of the initial guesses for the parameters were 10m4, 10m9 and 10-14, with initial step sizes of lo-“, lo-l5 and lo-*’ respectively. For the optimum escalation of the iteration two values were found, lo3 and lo’, which both gave the minimum value 0.12154 for the function F; the lower value (103) was chosen for the final calculations. The optimization of the system suggests that the error function is very complicated, so the observed data were treated by graphical analysis.” The error function used in the iteration process is drawn as contour plots and three-dimensional surfaces in Figs. 4 and 5. The function F splits into many minima, but the global minimum is clearly distinguishable from the local ones. The optimization of the iteration is necessary, otherwise the iteration may lead into a wrong pit and give erroneous results. In Table 4 the results calculated by the three different methods are collected together. The calculated p/?- and pK-values are very close to each other and stay within the error limits. In the matrix method the error of the first dissociation constant is smaller than in the other methods, but the errors of the second and third constants are one and a half times
6
Fig. 2. “C NMR
curve,
of the acid CISHz60.+
464
ANNE-MARIE SALONEN
Table
1. Effect of the initial guesses of the dissociation constants when respecttvely for each parameter and the escalation Initial
guess
PI
82
10-j ]O-4 10-z lo-6 lO-4 lo-4 lO-4
IO-” IO-’ 10-‘O lo-” 5 x 10m9 4 X lo-’ 5 X lOmY
Computed 104 x 8,
P? ._ IO-” lom’J lo-l5 10-lh lo-” ; ;;I:: E
lo9 x Bz
1 652 1 255 1 250 1.130 1.219 1.652 1.614
the inittal step size is IO-“‘, of the iteration is lo3
lOI x &
pK,
PK,
PK,
F
3.069 2.360 2.361 2.143 2.301 3.069 2.997
3.782 3.901 3.903 3.947 3.914 3.782 3.792
4.525 4 517 4 517 4.514 4.516 4.524 4.524
5.207 5.208 5.208 5.208 5.208 5.207 5.207
0.47827 0.12154* 0.12247 0.16241 0.12586 0.47827 0.41325
4.940 3.812 3.810 3.450 3.712 4.940 4.828
step size
Computed
PI
BZ
B,
IO4 x B,
10’ X 81
10’” x &
pK,
PK,
lo-‘0 lo-” 10-I’ 10-P lo-l4 lo-‘5 10-I” IO-” lo-‘* lom’v
,0-u lo-‘0 10-J’ 10-18 ;;I:::
0.5072 0.5074 0.5074 0.5087 1.239 1.255 1.227 1.249 1.178 1.000
1.693 1.693 1.693 1.693 3.767 3.812 3.73 I 3 795 3.593 1.452
1.041 1.041 1.041 1.041 2.332 2.360 2 310 2.349 2.224
4.295 4.295 4.295 4.293 3.907 3.901 3.911 3.903 3.929 4.000
4.476 4.477 4.477 4.478 4.517 4.517 4.517 4.517 4.515 4.838
10-24
for each parameter
parameters
lo-5 lo-6 IO-’ 10-R lo-9 lo-‘0 lo-” ,0-l> lo-” lo-‘4
]O_” ,0-X 10-P
1O-*0
parameters
Table 2. Effect of the inittal step size when Initial guesses are 10m4. 10m9 and lo-l4 respectively and the escalation of the Iteration is 10’ Imttal
lO-‘5 and
1.ooo
F ..~ 1.46573 1.46567 1.46565 1.46544 0.12254 0.12154* 0.12421 0.12178 0.13777 44 22171
pK3 5.211 5.21 I 5.211 5.211 5.208 5.208 5.208 5 208 5.208 5.162
Table 3. Effect of the escalation of the iteration when imtial guesses are 10m4. lO-9 and lo-“’ and initial step stzes are IO- “’. lO_” and lo-” respectively for each parameter Computed Escalatton 10 10’ 10’ IO4 IO5 lo6 10’ 10X 10Y 10”’
parameters
IO4 x /I,
IO9 X /I2
lOI x /I?
PK,
PK,
PK,
F
1.249 1.227 1.255 1.239 1.255 1.087 1.224 1.075 1.177 1 149
3.795 3.731 3.812 3.767 3.810 3.335 3 722 3.299 3.590 3.510
2.349 2.310 2.360 2.332 2.359 2.063 2.304 2.041 2.222 2.172
3.903 3.91 I 3.901 3.907 3.901 3.964 3.912 3.969 3.929 3.940
4.517 4.516 4.517 4.517 4.517 4.513 4.517 4.513 4.516 4.515
5.208 5.208 5.208 5.208 5.208 5.209 5.208 5.209 5.208 5.208
0.12178 0.12421 0.12154* 0.12254 0.12154* 0.19335 0.12477 0.20403 0.13827 0.15130
log F
Fig. 4. Contour plot of the logarithm of the error function as a function of the cumulative drssociation constants /I, and 8,: fi3 IS fixed at the opttmal value 2.36 x 10m’4.
Fig. 5. Three-dtmenstonal surface for the logarithm of the error function. The lowest contour hne is -0.95 and the htghest 5.35: the distance between the lines IS 0.2.
Table 4. Cumulative
Dissociation
constants
of a tribasic
organic
acid
and stepwise
dissociation constants, with standard with different programs
465 devtattons
computed
Programs Parameters ;, ; s;;;; sl f &,)
KDISS 1.25, x 1o-4 3.82, x 1O-9 2.36, x IO-l4
Matrix
method
(1.25,+0.11)X 10-a (3.82, _+ 0.32) x 1O-9 (2.36, _+ 0.20) x IO-l4
PB, + s(pP,) P/lz + s(p82) P& * S(PP,)
3.901 * 0.041 8.418 + 0.018 13.627 + 0.013
K, f s(K,)
1.25, x IO-” 3.03, X 10-S 6.19, x 10-h
(1.25,f_O.17) x IO-” (3.03, +_ 0.16) x IO-’ (6.19, & 0.05) x lo-”
3.03, x 1om5 6.19, x 10m6
3.901 4.517 5.208
3.900 * 0.059 4 517 * 0.022 5.208 +_ 0.035
3.899 4.517 5.208
Kz + s(Kz) Kj + s(K,) PK, * ~(PK,) pK, * .s(p&) PK, +_s(pK,)
as great that the mization. iteration inverting a large requires
NDISS (1.25,f0.12) X 10-d (3.81? 5 0.18) x 1O-9 (2.36,+ 0.13) x lo-l4
as in the iteration
method.
It is thus
3.900 8.417 13.626
obvious
iteration method is quite useful after optiThe matrix method is more rapid than the method, but requires a computer capable of a matrix with its elements distributed over number of orders of magnitude, and also great precision of calculation.
Acknow/edgemenr-The and Miqa Samppala
author IS grateful to Timo Nurmi for recording the “C NMR spectra. REFERENCES
1. A. Albert and E. P. Serjeant, The Determination of Ionization Constants, 3rd Ed., p. 61. Chapman and Hall. Cambridge, 1984. 2. R. F. Cookson, Chem. Rev.. 1974, 74, 5.
3.899 + 0.039 8.417 + 0.036 13.625 _+ 0.037
1.25,x 10-a
A. Loewenstein and S. Meiboom, 3. E. Grundwald, J Chem. Phys.. 1957. 27, 641. 4. J. F Hinton and E. S. Amis, Chem. Rev., 1967,67,367. Minimum of a Function of Several 5. J. P. Chandler, C’ariables. Program 66.1. Quantum Chemistry Program Exchange, Indiana University, Bloomington, Ind., 1966. Blometrika, 1956, 43, 353. 6. M. H. Quenouille, Multivariate Error Analysis. p. 30. 7. A. A. Clifford, Apphed Science Publishers, London, 1973. 8. Y. Bard. Nonlinear Parameter Estimation, p. 290. Academic Press, New York. 1974. 9. N. Draper and H. Smith, Applied Regression Analysis, 2nd Ed., p. 70. Wiley, New York, 1981. 10. A. Loewenstein and J. D. Roberts, J. Am. Chem. Sot., 1960, 82, 2705 University of Turku. Unpublished pro1 I. S. Pihlajamaki. grams DACONT and THREE/Z.