Microcomputer application of non-linear regression analysis to metal-ligand equilibria

7’0hraVol. 35, No. 7, pp. 507412, 1988 Rintcd in Great Britain. All rights reserved

Copyright 0

0039.9140/88 S3.00 + 0.00 1988 Pcrgamon Pma plc

MICROCOMPUTER APPLICATION OF NON-LINEAR REGRESSION ANALYSIS TO METAL-LIGAND EQUILIBRIA PAULD. TAYLOR*,IAN E. G. MORRISON and ROBERT C. HIDER Department of Chemistry and Biological Chemistry, University of Essex, Wivenhoe Park, Colchester, England (Received 5 October 1987. Revised 14 January 1988. Accepted 17 February 1988)

Summary-A non-linear least-squares regression program is described which is suitable for PC-compatible microcomputers. The program is written in GWBASIC, but compiled to run with the Intel 8087 fast numeric processor. Subroutines which simulate functions are compiled separately from the main program. Parameters are optimized by a Gauss-Newton-Marquardt algorithm which can be provided with either analytically or numerically calculated partial derivatives. Multi-component potentiometric titrations are simulated and parameters optimized by using analytical derivatives. Spectrophotometric titrations are also simulated, but absorptivities are optimized by linear regression while stability constants are optimized non-linearly by using numerical derivatives. Provision is made for “global analysis” of parameters. The experimental points can be displayed on screen, along with the “best” fit and the speciation. The program is demonstrated here by the determination of the pK, values and stability constants of a hydroxypyridinone ligand and its complexes with Fe(III).

Despite an established library of powerful mainframe computer programs for the determination of stability constants from titrimetric data,’ efforts have been made by numerous workers” to optimize stability constants by use of microcomputers. The impetus for this is convenience, especially with the increasing use of microcomputers for the control of automated systems, in which data are saved on disc and directly input to data-processing programs in the laboratory. Despite its faults, BASIC is a good language for such programs, because it is widely understood and accessible to adaptation by those who are not computer experts. The computer program NONLIN15, written in Microsoft GWBASIC, was conceived as a laboratory tool for the determination of stability constants of multi-component equilibria from data gathered by an automated titrator. This system, which will be described in a later paper, simultaneously obtains equilibrium readings of glass electrode emf and optical absorbance at a single wavelength, for incremental titrant additions. To speed up NONLINIS it is compiled to operate on a PC-compatible computer (Opus PCII) equipped with an Intel 8087 fast numeric by the Microway 8087 BASIC/ processor, INLINE software Version 1.15, using its facility for separately compiled subroutines, which allows new functions to be included with minimal effort. A speed that is more than SO times that with interpreted BASIC is achieved.

PROGRAM STRUCTURR Simulation of functions

The separately compiled subroutines which simulate chosen functions perform the input of all function parameters, calculation of statistical weights, generation of y as a function of x, and calculation of analytical partial derivatives of this function with respect to specific parameters. Simulation of potentiometric titrations The chemical system modelled consists of a metal ion and ligand in equilibrium with a series of complexes which tend to dissociate as the hydrogen-ion concentration is increased. The computer simulation of this model involves the solution of the set of simultaneous equations representing the conservation equations and equilibrium quotients. The free metal ion, ligand and hydrogen-ion concentrations in this system are calculated by Newton-Raphson iteration. The three appropriate root equations shown below use the following terms:

B(j) S(j) 4&j) Ct(i)

c(i) *Author for correspondence. 507

Cumulative formation constant for the jth associated species. Concentration of the jth associated species. Stoichiometric coe.ti%ient of the ith component in the jth associated species. Total analytical concentration of the ith component. Free concentration of the ith component.

508

PAUL

Mt, Lt, Ht

m,4h NC NS Kw 60’)

D.

TAYLOR

Total analytical concentrations of metal ion, ligand and hydrogen ion respectively. Free metal ion, ligaud and hydrogen-ion ~n~ntrations respectively. Total number of components. Total number of associated species. Ionization constant of water [H+J[OH-1, Absorptivity of the jth species NS

1

Mt = m f

e(m, P(j)

J-1

(1)

Lt = I + Ce(L iMY)

(2)

Ht =ih + xe(h, j) S(j) - Kw/h

(3)

et at.

than 0.01% of the analytical concentrations Mt, Lt and Ht. At this point the concentrations of all the complex species and the value of -log h corresponding to a given volume of acid or base titrant will have been calculated from the current estimates of the stability constants. Also at this point, the Jacobian matrix is fully inverted by using the last steps of the Choleski algorithm. The inverse Jacobian thus generated is used, as described by Avdeef and Raymond,‘O to generate. the analytical partial derivatives required for the Gauss-Newton-Marq~rdt optimization algorithm. Values for 8 log h/a~t, a log h/dLt and 8 log h/aHt are obtained directly from the inverse Jacobian, whereas terms for a log h/a log B(j) are obtained from equation (9).

where S(j) = p(j) *$ c(i)d’*fi

(4)

The 3 x 3 Jacobian matrix of partial derivatives9 of each total analytical concentration Mt, Lt and Ht with respect to the logarithm of the concentration of the free components m, I and h [equations (5) and (6)] has been shown to be a symmetrical matrix,‘0 so calcuIation is reduced by one third. The use of logarithmic space also helps to avoid computer overtlow and underflow errors.

The inverse Jacobian matrix is also used to extrapolate the titration curve to give better initial estimates for free component concentrations at the next point.” The statistical weighting scheme suggested by Avdeef” is used. Rigorous statistical weighting is vital when the sum-of-squares of emf or pH is minimized. Simulation of spectrophotometric titrations Calculation of absorbance as a function of -log h involves solution of one less simultaneous equation than for potentiometric titrations. A 2 x 2 Jacobian matrix is therefore employed. In other respects the algorithm to generate speciation is identical to that already described. Absorbance is then calculated from equation (10).

P

m + f

e(m, j)%(j)

J=l

J=

~4mAeGjPO’)

1

etm,j~{~,j)S{j)

‘

The vector X~+ 1giving the next improved estimate of log m, log I and log h is calculated from the previous estimate by equations (7) or (8) where J is the Jacobian matrix and F, is the column vector consisting of the current solution of the root equations.

J&+I -X,,)=

-F,

%+1 =K-

(7) J-IF;,

(8)

Solution of ~mul~neous equations by full matrix inversion as shown in equation (8) is computationally inefficient9 and unnecessary at this stage. Choleski’s method9 is employed instead. Iteration continues until the equations are solved with an accuracy better

(10) The analytical derivatives for absorptivities supplied to the Gauss-Newton optimization algorithm are simply the relevant species concentrations as shown in equation (11). aA

-=S(j)

WI

The remaining parameters describing spectrophotometric titrations are optimized by using numerically obtained derivatives for simplicity (see next section).

509

Microcomputer non-linear regression analysis Optimization of parameters by the Gauss-NewtonMarquardt algorithm Minimization of the sum-of-squares of residuals is performed by a standard Gauss-Newton algorithm developed from the program GENFIT.” The Marquardt algorithm14 has recently been introduced into the program and functions by multiplying the diagonal terms of the Hessian matrix by the factor (1 + x). The value of x is initially 0.1 and is reduced by a factor of three each time the sum-of-squares of residuals is reduced in an iteration. When convergence fails, the parameters of the previous iteration are restored and x is repeatedly increased by a factor of three until convergence is regained. Up to ten variables may be refined simultaneously. Some may be held constant for a sequence of operations and freed for variation later. If analytical derivatives have not been supplied by the function simulation algorithm than a standard numerical differentiation method is employed. Once convergence has been established (tested by the ratio of successive iterations) the statistical parameters chi-squared and the Hamilton R factori are calculated. In both potentiometric and spectrophotometric analyses the stepwise stability constants are optimized rather than the overall stability constants. A special modification to the Gauss-Newton algorithm has been made to deal with linear parameters such as absorptivities. The method is based on that reported for the program ELORMA.3 Numerical derivatives for non-linear parameters such as stability constants are calculated by the usual perturbation method, but immediately after the speciation has been generated the absorptivities are replaced by their best linear least-squares values. This separate treatment of linear and non-linear parameters gives more reliable and more rapid convergence.4 Zuberbiihler et al. have employed a highly sophisticated mathematical technique in their program SPECFIT which enables the use of analytical derivatives for non-linear parameters to be retained4 during this two-stage process. Global optimization Global optimization of the parameters common to a set of experiments gives more reliable values and a better chance of rejection of a spurious model function. Parameters such as equilibrium constants must by definition (assuming ideality) be constant and therefore common to a set of titrations. For parameters which are not constant throughout a set of titrations, another method of global optimization is required. The method used recognizes the experimental ease with which some variables such as concentration may be increased by a known factor even when the absolute values are not known. For example, in a set of titrations in which the metal concentration is increased between titrations by addition of 1,2,3 ml, etc., of a stock solution, the factors

1, 2, 3, etc., may be input into NONLINlS and held constant while the metal concentration is optimized. Experimentally this is performed by repetitive titration of a single sample solution but with addition of say 1 ml of stock metal or ligand solution between titrations. Dilution effects from these additions and from the titrant volume are all handled by the program. This process can help to eliminate small systematic errors. Utilization of combined potentiometric photometric data

and spectro-

Programs which perform non-linear regression analysis simultaneously on potentiometric and spectrophotometric data have been described in the literature,‘” but this approach requires the weighting of the two types of data to be in the correct proportion. A different approach is used in NONLINIS, in that the two sets of data are analysed in a two-stage process which exploits the properties of each. The parameters for initial acid concentration, titrant concentration, pKw and carbonate concentration have no effect on the plot of absorbance vs. Nemstian electrode response. Furthermore, the process of resolving the absorbance into the contributions made by each absorbing species is independent of the electrode E,,. The spectrophotometric data may therefore be used in isolation to optimize speciation. An arbitrary E0 value may then be used to give apparent but erroneous values of the equilibrium constants. The way in which errors in E,, are transformed into errors in these apparent equilibrium constants is shown in equation (14). This function is dependent on the proton stoichiometry as well as the error E, in the arbitrary E,,. The automated spectrophotometric/potentiometric titrator described earlier gathers the data for emf vs. titrant volume point by point, simultaneously with the spectrophotometric data. The potentiometric data are then used to perform an internal calibration by optimization of E, and adjustment of the equilibrium constants according to equation (14), while the absorbing species concentrations are held constant at the levels determined spectrophotometrically. The initial guesses for this second stage are the arbitrary E, used in the analysis of the spectrophotometric dam and the corresponding apparent equilibrium constants. The analytical derivative supplied to the Gauss-Newton-Marquardt algorithm is shown in equation (16). log S(j) = log /3(j) + e(m, j) log m + e(l, j) log 1 + e(h, j) log h emf=E,,+E,+Slogh

(12) (13)

where S is the electrode slope. For fixed speciation: logfl(j)

= const + e(h, j)(emf - E0 - E,)/S d log /?(j)/d E, = - e(h, j)/S

(14) (15)

PAUL D. TAYLOR et al.

510

Table 1. Comparison of data” analysed by MINIQUAD and NONLINI 5; the standard error of the emf was taken as 0.6 mV, that of the titrant volume as 0.010 ml; these values were substituted into the weighting equation of Avdeef12 Titration

log BOJ,,

MINIQUAD~

Program

5 6 7 8 9

9.659

NONLINIS

5 6 7 8 9

log BW 12.071

‘og &,r 5.636 (8)* 5.646 (14) 5.614 (15) 5.635 (3) 5.616 (2)

log 81.02 10.396 (7)* 10.408 (14) 10.368 (14) 10.391(4) 10.366 (6)

log 81.0.X 13.824(13)* 13.849 (21) 13.736 (22) 13.811 (25) -

5.636 (11) 5.648 (11) 5.614 (10) 5.634 (2) 5.615 (2)

10.396(11) 10.407(11) 10.367 (9) 10.392 (4) 10.361(10)

13.822(15) 13.848(10) 13.745 (9) 13.797 (2) -

*The digits in parentheses are IO3times the standard deviations of the mantissas of the log /? values.

From lotion d emf -= dE,

(13) 1+&S? I=I

alogh

a himi).

standard error.” Convergent for titration 7, consisting of 97 points, occurred from initial estimates of log B,,O,,= 5, log /?1,0,2 = 9 and log /?,,O,S = 12 in 3 iterations, taking a total time of 9 sec. AppIication

The partial derivative 8 log h/a log flu) can be found from the inverse Jacobian matrix as already described. The remaining parameters for initial acidity, titrant concentration, pKw and carbonate concentration may be optimized along with E, if desired. A model must give an acceptable fit with both the s~trophotomet~c and potentiomet~c data if it is not to be rejected. The method is found to be particularly advantageous for measuring dissociation occurring at pH values greater than 10 or less than 4. In these circumstances potentiometric titration becomes less reliable as 8 log h/a log @G) g 1. The method afso avoids independent optimization of highly correlated variables such as initial acidity and log flu) for a minor species undergoing dissociation at pH values near 7. The technique can be used to optimize E, in titration of meal-li~nd systems. However, it is not then possible to hold constant any plya values which produce significant speciation change over the titration range. The method demands that all appreciable proton-dependent speciation be optimized simultaneously from the spectrophotometric data. Validation

NONLINIS has been validated by repeating the analysis of the potentiometric titration data of the nickel-glycinate system studied by Braibanti et al.” NONLINlS optimizes the logarithms of stepwise stability constants, which are summed in Table 1 for comparison with the overall logarithms of the stability constants obtained with MINIQUAD. The standard errors given here for NONLINlS are similarly obtained by summation of the stepwise standard errors. This gives a pessimistic estimate of the true

The use of the program is briefly demonstrated here for determination of the stability constants of the bidentate ligand 3-hydroxy-2-methyl-4” (Ill)-pyridinone (I). 0

I

Nitric acid (0.2M) and potassium hydroxide solution (0.2M) were prepared from BDH “Convol” ampoules. Ionic strength was maintained at 0.2 with potassium nitrate (BDH Aristar grade). All solutions were made up with 18-MR/cm water from a Millipore Milli-Q system. Temperature was maintained at 22.5 f 0.1” with a Churchill t~e~~~~ator. Ambient temperature was also kept at 22.50. The ferric ion solution was prepared from the Fisons analytical grade nitrate and standardized by EDTA titration according to established procedures.‘9 The ligand 3-hydroxy-2-methyl-4( lH)-pyridinone was synthesized and purified in this department.” Glass electrode performance was first checked by titrating with potassium hydroxide solution (0.2M) a 25.ml volume of potassium nitrate solution (0.2M) acidified with nitric acid, and making a Grans piot.2’ The &and solution was then added, and the plya values were determined from three titrations, without removal of the electrode, each performed after addition of a further 2 ml of nominally 5mAf ligand somtion. Absorbance was monitored at 270 nm. In a separate experiment 5 ml of nominally 5mM &and and 1 ml of nominally 5mM ferric nitrate were added to 25 ml of potassium nitrate sohttion (0.2&f) before titration with potassium hydroxide (0.2M). Subsequent titrations were performed after addition of a further 1 ml of ferric nitrate solution, until a total of 3 ml had been added. Absorbance was monitored at a wavelength of 305 mn, where the uncomplexed ligand has a low absorbance. compared with the complex.

Microcomputer

4.2

1.13 x lo-‘,

I

I

I

49

5.6

6.3

-emf /59 16

-emf/59.16

Fig. 1. Spectrophotometric titrations of the ligand (I). Absorbance at 270 nm us. measured -emf/59.16. A further 2 ml of ligand solution was added between titrations. The solid line is the function simulated from the globally optimized parameters below. E, was fixed at zero. Iteration 5, Chi-squared = Error limit 1.1116 x 10m3 f for LH = 6 for LH: = Absorbance zero = pKa =

511

non-linear regression analysis

R = 1.16 x 10m3

Fig. 3. Spectrophotometric titration of ligand (I) in the presence. of 1, 2 and 3 ml of ferric nitrate solution. Mt : Lt = (1) 1: 5.809; (2) 1: 2.905; (3) 1: 1.936. The solid line is the function simulated from the globally optimized parameters below. Iteration 5, Chi-squared = 7.59 x lo-‘, R = 1.36 x IO-” Error limit 7.3509 x IO-’

1262.8 f 0.5 l.mole-‘.cm-’ 796.3 f 0.5 I.mole-‘.cn-’ 0.3119 f 0.003 4.609 f 0.001

L for ML2+ = L for ML: = 6 for ML, = t for LH = 6 for LH: = Log K2 = Log K3 =

RESULTS

Experimental points and the simulated curves generated by the computer with optimized parameters are shown in Figs. 1-3, and the speciation at convergence is shown in Fig. 4. The optimized parameters are given in the legend to each figure. The slope of the glass electrode response as a fraction of the theoretical Nemstian response at 25” was optimized from spectrophotometric titration of a highly purified ref-

290+ 1 l.mole-‘. cm-‘. 674.7 f 0.6 l.mole-I. cm-‘. 1354& 1.4 l.mole-‘. cm-‘. 118.2 f 0.8 l.mole-I. cm-‘. 93.2 + 0.5 I.mole-‘. cm-‘. 12.027 rf:0.001 9.734 f 0.001

erence ligand with a single pKa (maltol). The pKa of the ligand 3-hydroxy-2-methyl-4( 1H)-pyridinone was first globally optimized from the spectrophotometric titration curves by using an assumed electrode E0 of zero. The error E,, in E, was then optimized from the simultaneously measured data for emf us. titrant volume [equations (14) and (1611. The stability constants for the complexes of the ligand with Fe(II1) were optimized from the spectrophotometric titration of the metal-ligand system, by using the pKa values, electrode slope and &value determined above. Comparison with other programs The program NONLIN15 has many similarities to the BASIC program MICMAC.* Both are in modu100%

4

Microlitres

of base

Fig. 2. Potentiometric titrations obtained simultaneously with the absorbance data shown in Fig. 1. The error ,T, in E, used for Fig. 1 was optimized by using the derivative shown in equation (16), and the pKa was adjusted according to equation (14). The solid line is the function simulated from the globally optimized parameters below. Iteration 8, Chi-squared = 0.370, R = 9.58 x lo-’ Error limit 0.3621 Initial acid cont. = 4.348 x lo-’ f 3 x lo-‘A4 E, = 0.9086 f 0.0005 Added acid=419fO.l ~1 pKa = 3.700 f 0.0015

4.2

4.9

5.6

6.3

-emf/59.16 Fig. 4. Speciation plot for the titration shown in Fig. 3, curve 2, in which 2 ml of ferric ntrate were added. Mt:L1=

1:2.905.

512

PAUL D. TAYLORet al.

lar form to simplify the inclusion of new model functions, but NONLIN15 can use analytical derivatives, and this speeds operation and improves the reliability of convergence. When potentiometric titrations are modelled these analytical derivatives are used to optimize the total component concentrations, &values and other parameters, as well as the logarithms of all the equilibrium constants. In this respect it is similar to the program TITFIT.* However, potentiometric parameters are optimized in NONLINIS by ~ni~ng the squared residuals of emf rather than titrant volume. A rigorous slopedependent weighting scheme is used. Minimization with respect to titrant volume inherently gives the equivalence regions less weight by virtue of the slope-dependent density of points, and this method or the method adopted in NONLINl5 have been found by other workers to give closely comparable statistical parameters.= A recent publication23 describes the FORTRAN program PROTAF with the features above, which minimizes the sum of square residuals of titrant volume and emf or pH. The presence of carbonate in the basic titrant is included in the ~tentiomet~ model used for NONLINl5. When NONLINlS is used to analyse spectrophotometric data, the optimization of absorptivities proceeds by a method similar to that used in the program ELORMA: as discussed previously. A recent pub~~tionz4 describes the FORTRAN 77 program CFTSP which optimizes stability constants and absorptivities, but without the distinction between the linear and non-linear parameters. Finally, a simple method of utilizing simultaneously measured spectrophotometric and potentiometric data has been included in the program. Satisfactory speed of operation on a ~crocomputer is achieved by use of the Intel 8087 coprocessor. Several improvements to the program are currently being considered. One of these is the inclusion of a general speciation program called MLGEN19,25 which is routinely used in this department to simulate speciation in any threecomponent system with up to thirty complexes. Currently, however, NONLINl5 offers a menu of fixed but comprehensive equilibrium models in which nonexistent species may be suppressed by use of large negative log K or pKa values.

The GWBASIC request, Acknowledgement-This

program

listing is available

research was supported

on

by the

British Technology Group.

REFERENCES D. J. Leggett (ed.), Computational Method for the Determination of Stability Constants, Chap. 1. Plenum, New York, 198.5,and references therein. 2, A. D. Zuberbiihler and A. T. Kaden, Talanta, 1982,29, 201. 3. H. Gampp, M. Maeder and A. D. Zuberbiihler, ibid., 1980, 27, 1037. 4. H. Gampp, M. Maeder, C. J. Mayer and A. D. Zuberbiihler, ibid.. 1985. 32. 257. 5. R. J. Motekaitis and At E. Martell, Can. J. Chem., 1982, 60. 168. 6. A.‘Sabatini, A. Vacca and P. Cans, Talanta, 1974,21, 53 (original mainframe proaram): A. P. Arnold. S. A. Daignault and D. L. Raben~ein,‘,&al. C/rem., 1985,57, 1112 (Microcomputer version). 7. A. Izquierdo and J. L. Bettran, Anal. C&n. Acta, 1986, 181, 87. 8. A. Laouenan and E. Suet, Tafanta, 1985,32, 245. 9. G. H. Golub and F. C. Van Loan, Matrix Computations, North Oxford Academic, Oxford, 1983. 10. A. Avdeef and K. N. Raymond, Inorg. Chem., 1979,18, 1605. 11. D. J. Leggett, op. cit., p. 295. 12. A. Avdeef, Anal. Chim; Acta, 1983, 148, 237. 13. I. E. G. Morrison, unnublished nroaram. 14. D. W. Marquardt; J. $0~. Znd. >ppL Math., 1963, 11, 431; W. H. Press, B. P. Flannery, S. A. Teukolsky and W. T. Vetterling, Numerica/ Recipes: The Art of Scientz$c Computing, Cambridge Uuiversity Press, Cambridge, 1986. 15. M. Jaskolski and L. Lomozik, Talanta, 1985, 32, 511. 16. D. J. Leggett, op. cit., p. 293. 17. A. Braibanti, F. Dallavalle and G. Mori, Ann. Chim., Rome, 1978, 68, 861. 18. M. Micheloni, A. Sabatina and A. Vacca, Inorg. Chim. Acfa, 1977, 25, 41. 19. A. I. Vogel, A Textbook of Quantitative Inorganic Analysis, 3rd Ed., Longmans, London, 1961. 20. 3-Hydroxy-2-methyl-4-(lH)-pyridinone was synthesized by P. Sarpong at the University of Essex. 21. G. Grans, Acta Chem. &and., 1950, 4, 559; Analyst, 1952, 77,661. 22. A. D. Zuberbiihler, Personal communication. 23. R. Fournaise and C. Pettifaux, Tafanta, 1987, 34,385. 24. L. Lampugnani, L. Meites, P. Papoff and T. Rotunno, Anal. Chim. Acta, 1987, 77, 194. 25. P. D. Taylor, Ph.D. The&s, University of Essex, 1987.