Determination of Equilibrium Constants with the OPIUM Computer Program

\underline{Mojmír Kývala}^{a}, Přemek Lubal^{b}, Ivan Lukeš^{a}

aDepartment of Inorganic Chemistry, Charles University, Hlavova 2030, 128 43 Praha 2, Czech Republic bDepartment of Analytical Chemistry, Masaryk University, Kotlářská 2, 611 37 Brno, Czech Republic

The recently published^{1} OPIUM computer program for the analysis of solution equilibria by the generalized least squares method is now compared with several well-known programs for the analysis of potentiometric, spectrophotometric (or enthalpimetric) and NMR chemical shift (or capillary zone electrophoresis effective mobility) data. To increase the objectiveness of the comparison, only “third party” data are treated (formation of polynuclear complexes, ligand competition, dissociation of a weak acid).

The main feature of the OPIUM program (besides extensive analysis of both the estimates and the residuals¶) is a strict separation of the dependence of a response variable y on the calibration and chemical parameters y = f(ϑ, c, n, V, c’, q, β, x) into two independent parts. These are a calibration function (the dependence of the response variable y on the calibration parameters ϑ and the free concentrations c of all components) y = ϕ(ϑ, c) and the system of nonlinear equations (mass balances)

c_{k} + ∑_{i=1}^{m} q_{ik}β_{i} ∏_{j=1}^{n} c_{j}^{q_{ij}} = \frac{c_{k}V + c’_{k}x}{V + x} = \frac{n_{k} + c’_{k}x}{V + x} , k = 1, …, n ,

n is the number of reactants, m is the number of complexes, c_{j} is the free concentration of the j-th reactant, q_{ij} is the element of the matrix of stoichiometric coefficients q ∼ (m, n), β_{i} is the stability constant of the i-th complex, c_{j} (n_{j}) is the analytical concentration (mass amount) of the j-th reactant in the titrand, V is the volume of the titrand, c’_{j} is the analytical concentration of the j-th reactant in the titrant and x is the added volume of the titrant (only if a titration is applied). This separation (together with a correct definition of the variance-covariance matrix V of the random vector of the dependent variable y) is a key to the treatment of many different types of data in any combination.

Formally, \underline{any} subset of the parameters ϑ, c, n, V, c’, q, β can be “marked” unknown and calculated. Moreover, simple boundary conditions for selected unknown parameters like for instance a_{1} = a_{2}+a_{3}† or c_{1} = 0.0736+4c_{2} ∧ c_{3} = c_{2} can be employed to suppress the correlation of the unknown parameters estimates. The OPIUM program also analyzes data containing sections with zero free concentrations of some reactants (necessary, for example, when calculating dissociation constants of a weak acid and stability constants of complexes of the acid with several metals together, i.e. from both types of titrations).

¶ including their bias and implied parametric and internal curvature of the model † a_{1}, a_{2}, a_{3} could be e.g. molar absorption coefficients

It is proved, that the OPIUM program is robust enough to mostly find the global minimum using automatically generated initial estimates of the unknown stability constants¶ (i.e. without any human assistance). This feature allows us to test \underline{many thousands} of chemical model hypotheses by the “trial and error” method in several hours† (the program is written in the standard ANSI C language and can be run on virtually any computer). The best model(s) can be, then, verified using the built-in ESI (\underline{e}stimation of \underline{s}toichiometric \underline{i}ndices) approach^{2}, i.e. a joint calculation of all the unknown stability constants and selected elements q_{ij} ∈ R of the matrix of stoichiometric coefficients q.

1 Description of the OPIUM Computer Program

1.1 Chemical model (mass balances)

To be a little more accurate, we have to define the number of reactants with unknown free concentration r ≤ n (a lot of spectroscopic data are collected using buffered solutions, i.e. the free concentration instead of the analytical concentration of H^{+} is known!). Then the mass balances are§

f_{k} = c_{k} + ∑_{i=1}^{m} q_{ik}c_{n+i} − C_{k} = 0 , k = 1, …, r ,

c_{n+i} = β_{i} ∏_{j=1}^{n} c_{j}^{q_{ij}} = β_{i}e^{∑_{j=1}^{n} q_{ij} ln c_{j}} = β_{i}10^{∑_{j=1}^{n} q_{ij} log c_{j}} = 10^{log β_{i}+∑_{j=1}^{n} q_{ij} log c_{j}} , i = 1, …, m .

Partial derivatives of the components of the of the function f with respect to the logarithms of the unknown free concentrations of reactants are

\frac{∂f_{k}}{∂ ln c_{k}} = 1 + ∑_{i=1}^{m} q_{ik}^{2}c_{n+i} , k = 1, …, r ,

\frac{∂f_{k}}{∂ ln c_{l}} = \frac{∂f_{l}}{∂ ln c_{k}} = ∑_{i=1}^{m} q_{ik}q_{il}c_{n+i} , k ≠ l = 1, …, r .

¶ If initial estimates of the unknown stability constants are not given, the OPIUM program generates them according to the simple formula log β_{i} = A + B ∑_{j=1}^{n} q_{ij} , A and B are empirical constants. If \underline{all} the other chemical parameters c, n, V, c’, q are known, the unknown parameters optimization method used usually ensures rapid convergence to the correct solution. If they are not, some other local minima often exist and better initial estimates of β might be required to find the “correct” one.

† while being somewhere else and/or doing something more exciting § let the n + m components are ordered such that the first r of them are reactants with unknown free (known analytical) concentrations, the next n − r are reactants with known free (unknown analytical)

concentrations and the last m are complexes

Partial derivatives of the function y = f(ϑ, c, n, V, c’, q, β, x) with respect to the unknown chemical parameters p ≡ c, n, V, c’, q, β are¶

\frac{∂y}{∂p} = ∑_{i=1}^{n+m} \frac{∂y}{∂c_{i}} \frac{∂c_{i}}{∂p} .

Derivatives \frac{∂c_{i}}{∂p}, i = 1, …, r, are the solution of the system of linear equations

\frac{∂c_{k}}{∂p} + ∑_{l=1}^{n} \frac{∂c_{l}}{∂p} \frac{1}{c_{l}} ∑_{{i=1 \atop q_{il}≠0}}^{m} q_{ik}q_{il}c_{n+i} = A_{k} , k = 1, …, r ,

for p ≡ c_{v}, n_{v}, V, c’_{v}

A_{k} = \frac{∂C_{k}}{∂p} †,

for p ≡ β_{u}

A_{k} = −q_{uk} ∏_{j=1}^{n} c_{j}^{q_{uj}} ,

for p ≡ q_{uv} ∧ k ≠ v

A_{k} = −q_{uk} ln c_{v} c_{n+u} ,

for p ≡ q_{uk}

A_{k} = −(1 + q_{uk} ln c_{k})c_{n+u} .

Derivatives \frac{∂c_{i}}{∂p}, i = r+1, …, n, are zero§, derivatives \frac{∂c_{i}}{∂p}, i = n+1, …, n+ m (i.e. \frac{∂c_{n+i}}{∂p}, i = 1, …, m), are

\frac{∂c_{n+i}}{∂p} = B_{i} + ∑_{{k=1 \atop q_{ik}≠0}}^{n} \frac{∂c_{k}}{∂p} \frac{1}{c_{k}} q_{ik}c_{n+i} , i = 1, …, m ,

¶ excluding two exceptions (see the description of calibration functions of types 1 and −5) † it is clear that

C_{k} = \frac{c_{k}V + c’_{k}x}{V + x} = \frac{n_{k} + c’_{k}x}{V + x} ,

\frac{∂C_{k}}{∂c_{v}} = \frac{∂C_{k}}{∂n_{v}} = \frac{∂C_{k}}{∂c’_{v}} = 0 , v ≠ k = 1, …, r ,

\frac{∂C_{k}}{∂c_{k}} = \frac{V}{V + x} , \frac{∂C_{k}}{∂n_{k}} = \frac{1}{V + x} , \frac{∂C_{k}}{∂c’_{k}} = \frac{x}{V + x} ,

(\frac{∂C_{k}}{∂V} )_{dc_{k}=0} = \frac{(c_{k} − c’_{k})x}{(V + x)^{2}} , (\frac{∂C_{k}}{∂V} )_{dn_{k}=dc_{k}V+c_{k}dV=0} = −\frac{n_{k} + c’_{k}x}{(V + x)^{2}}

§ This is not true, if the “known” free concentration of the i-th reactant is considered an unknown parameter! The OPIUM program provides this option, if an \underline{alternative definition} of a solution is used (solutions are not prepared by a titration).

for p ≡ β_{i}

B_{i} = ∏_{j=1}^{n} c_{j}^{q_{ij}} ,

for p ≡ q_{ik}

B_{i} = ln c_{k} c_{n+i} ,

B_{i} = 0 for any other p. 1.2 Measurements (calibration functions) Overview of the supported calibration functions follows¶. Type 1

y = E_{0} + S log c_{1} + j_{1}c_{1} + j_{2} \frac{β_{1}}{c_{1}}

ϑ ≡ E_{0}, S, j_{1}, j_{2}

\frac{∂y}{∂c_{1}} = \frac{S}{ln10 · c_{1}} + j_{1} − j_{2} \frac{β_{1}}{c_{1}^{2}} , \frac{∂y}{∂β_{1}} = \frac{∂y}{∂c_{1}} \frac{∂c_{1}}{∂β_{1}} + \frac{j_{2}}{c_{1}}

Note: Function of type 1 describes the dependence of the difference y of potentials of a “hydrogen” (glass) electrode and a reference electrode in a cell with non-negligible liquid-junction potential on the free concentration c_{1} (β_{1} ≡ K is the ionic product of water). Thus, H^{+} must be the 1st reactant (1st column of the matrix q), OH− must be the 1st complex (1st row of the matrix q).

Type 2

y = E_{0} + S log c_{1}

ϑ ≡ E_{0}, S Type 3

y = E_{0} + S log c_{k} + j_{1}c_{1} + j_{2} \frac{β_{1}}{c_{1}}

ϑ ≡ E_{0}, S, j_{1}, j_{2}; k ∈ ⟨1, n⟩ Type 4

y = E_{0} + S log c_{k}

ϑ ≡ E_{0}, S; k ∈ ⟨1, n⟩ Type 5

y = E_{0} + S log \frac{c_{k}}{c_{l}} + j_{1}c_{1} + j_{2} \frac{β_{1}}{c_{1}}

ϑ ≡ E_{0}, S, j_{1}, j_{2}; k ≠ l ∈ ⟨1, n⟩ Type 6

y = E_{0} + S log \frac{c_{k}}{c_{l}}

¶ including some “non-trivial” derivatives

ϑ ≡ E_{0}, S; k ≠ l ∈ ⟨1, n⟩

\frac{∂y}{∂c_{k}} = \frac{S}{ln10 · c_{k}} , \frac{∂y}{∂c_{l}} = −\frac{S}{ln10 · c_{l}}

Note: Functions of types 5 and 6 describe the dependence of the difference y of potentials of a redox (platinum) electrode and a reference electrode on the ratio of free concentrations c_{k}, c_{l} of two reactants (i.e. the simplest case of a redox system, e.g. Fe^{3+}/Fe^{2+}, Ce^{4+}/Ce^{3+} etc.).

Type 7

y = ∑_{i=1}^{n+m} a_{i}c_{i} = a^{T}c

ϑ ≡ a_{i}, i = 1, …, n + m

\frac{∂y}{∂a_{i}} = c_{i} , \frac{∂y}{∂c_{i}} = a_{i} , i = 1, …, n + m

Type 8

y = d ∑_{i=1}^{n+m} a_{i}c_{i} = d a^{T}c

ϑ ≡ a_{i}, i = 1, …, n + m; d ∈ R

Note: Functions of types 7 and 8 describe the dependence of an \underline{extensive} linearly additive quantity y (absorbance in some “slow exchange” spectrometry (in the case of inseparable peaks), usually UV/VIS spectrophotometry, reaction heat in enthalpimetry, conductivity in conductimetry etc.) on the free concentrations c of all components.

Type 9

y = \frac{c_{i}}{c_{j}}

i ≠ j ∈ ⟨1, n + m⟩

\frac{∂y}{∂c_{i}} = \frac{1}{c_{j}} , \frac{∂y}{∂c_{j}} = −\frac{c_{i}}{c_{j}^{2}}

Type −1

y = − log c_{1}

Type −2

y = − log c_{k}

k ∈ ⟨1, n⟩

Type −3

y = αc_{i}

ϑ ≡ α; i ∈ ⟨1, n + m⟩

Type −4

y = dαc_{i}

ϑ ≡ α; d ∈ R, i ∈ ⟨1, n + m⟩

Type −5

y = ∑_{i=1}^{n+m} a_{i}p_{i}c_{i}/C = a^{T}Pc/C ¶†

P = diag(p_{1}, …, p_{n+m})

ϑ ≡ a_{i}, i = 1, …, n + m

\frac{∂y}{∂a_{i}} = p_{i}c_{i}/C , \frac{∂y}{∂c_{i}} = p_{i} (a_{i}C − ∑_{j=1}^{n+m} a_{j}p_{j}c_{j}) /C^{2} = p_{i} (a_{i} − y) /C

i = 1, …, n + m

\frac{∂y}{∂q_{uv}} = ∑_{i=1}^{n+m} \frac{∂y}{∂c_{i}} \frac{∂c_{i}}{∂q_{uv}} + p_{v}c_{n+u} (a_{n+u} − y) /C , u = 1, …, m , v = 1, …, n

Note: Function of type −5 describes the dependence of an \underline{intensive} linearly additive quantity y (chemical shift or relaxation rate, i.e. reciprocal relaxation time, in “fast exchange” NMR spectroscopy, effective mobility in capillary zone electrophoresis of a labile system etc.) on the free concentrations c of all components.

Type −6

y = p_{i}c_{i}/C ¶†

i ∈ ⟨1, n + m⟩

\frac{∂y}{∂c_{i}} = p_{i}(C − p_{i}c_{i})/C^{2} = p_{i}(1 − y)/C

\frac{∂y}{∂c_{j}} = −p_{j}p_{i}c_{i}/C^{2} = −p_{j}y/C , j = 1, …, n + m ∧ j ≠ i

\frac{∂y}{∂q_{uv}} = ∑_{j=1}^{n+m} \frac{∂y}{∂c_{j}} \frac{∂c_{j}}{∂q_{uv}} + A_{uv}

¶

p_{i} ∈ N , i = 1, …, n , p_{n+i} = ∑_{j=1}^{n} p_{j}q_{ij} , i = 1, …, m , C = ∑_{i=1}^{n+m} p_{i}c_{i} , \frac{∂C}{∂q_{uv}} = p_{v}c_{n+u}

† A little more complicated definition (using almost general diagonal matrix P) enables the correct description of the dependent variable no matter what linear transformation of the matrix q is actually in use. Moreover, dynamic definition of the analytical concentration C of the “active” component (not exclusively reactant or complex) is \underline{necessary} for the calculation of unknown chemical parameters c, n, V, c’, q (the value of C changes).

A_{uv} = p_{v}c_{i} (1 − y) /C , n + u = i , v = 1, …, n

A_{uv} = −p_{v}c_{n+u}y/C , u = 1, …, m ∧ n + u ≠ i , v = 1, …, n

Note: Function of type −6 describes the dependence of an \underline{intensive} quantity y (intensity of an isolated (separated) peak in some “slow exchange” spectrometry, usually NMR spectroscopy, intensity of an isolated peak in capillary zone electrophoresis of an inert system etc., normalized to the sum of intensities of all the peaks embracing \underline{all} components containing the active component) on the free concentrations c of all components (but mainly on c_{i}).

1.3 Regression model

Let m is the number of points, i.e. the number of ordered pairs of measured variables on abscissa and ordinate [x, y] ([x_{i}, y_{i}], i = 1, …, m) and n is the number of unknown parameters. If m = n, the values of the unknown parameters are the unique solution of a system of equations. If m > n, the system of equations is over-determined. Moreover, if the values of variables x, y and known parameters p ≡ ϑ, c, n, V, c’, q, β are accompanied with measuring errors, the solution of the system of equations is not unique.

Let the measuring errors ε_{x}, ε_{y}, ε_{p_{k}} of variables x, y, p_{k} are random, independent, have normal distributions, zero means, variances σ_{x}^{2}, σ_{y}^{2}, σ_{p_{k}}^{2} and act additively. Then the variables x, y, p_{k} are random, independent, have normal distributions, unknown means and variances σ_{x}^{2}, σ_{y}^{2}, σ_{p_{k}}^{2}

E(ε_{x}) = E(ε_{y}) = E(ε_{p_{k}}) = 0 ,

E(ε_{x}ε_{y}) = E(ε_{x})E(ε_{y}) = 0 ,
E(ε_{x}ε_{p_{k}}) = E(ε_{x})E(ε_{p_{k}}) = 0 ,
E(ε_{y}ε_{p_{k}}) = E(ε_{y})E(ε_{p_{k}}) = 0 ,
E(ε_{p_{k}}ε_{p_{l}}) = E(ε_{p_{k}})E(ε_{p_{l}}) = 0 , k ≠ l ,
E(ε_{x}ε_{y}) = E(ε_{x}ε_{p_{k}}) = E(ε_{y}ε_{p_{k}}) = E(ε_{p_{k}}ε_{p_{l}}) = 0 , k ≠ l .

However, it is usually assumed that only the values of the variables x and p_{k} are set (measured) independently in the course of a titration. Since it still holds that y = f(ϑ, c, n, V, c’, q, β, x) ≡ f(p, x), an error of the functional relation f(p, x) should be added to the “independent” measuring error ε_{y} of the variable y to form the true regression error ε.

From the linear Taylor approximation of the function f(p, x) in a close neighborhood of the means of the variables x and p_{k} it follows that

ε ≐ \frac{∂y}{∂x} ε_{x} + ε_{y} + ∑_{k} \frac{∂y}{∂p_{k}} ε_{p_{k}} ¶ ,

\frac{∂C_{k}}{∂x} = \frac{(c’_{k} − c_{k})V}{(V + x)^{2}} = \frac{c’_{k}V − n_{k}}{(V + x)^{2}}

¶ all the necessary derivatives are given, except for

E(ε) ≐ 0 , D(ε) ≐ V .

The elements of the variance-covariance matrix V are¶

v_{ii} ≐ D(ε_{i}) ≐ E(ε_{i}^{2}) ,

v_{ii} = (\frac{∂y}{∂x})_{x=x_{i},y=y_{i}}^{2} σ_{x_{i}}^{2} + σ_{y_{i}}^{2} + ∑_{k} (\frac{∂y}{∂p_{k}})_{y=y_{i},p_{k}=p_{k}_{i}}^{2} σ_{p_{k}_{i}}^{2} ,

v_{ij} ≐ cov(ε_{i}, ε_{j}) ≐ E(ε_{i}ε_{j}) ,
E(ε_{x_{i}}ε_{x_{j}}) = E(ε_{x_{i}})E(ε_{x_{j}}) = 0 , x_{i} ≢ x_{j} ,
E(ε_{x_{i}}ε_{x_{j}}) = σ_{x_{i}}^{2} , x_{i} ≡ x_{j} ,
E(ε_{y_{i}}ε_{y_{j}}) = E(ε_{y_{i}})E(ε_{y_{j}}) = 0 , i ≠ j ,
E(ε_{p_{k}_{i}}ε_{p_{k}_{j}}) = E(ε_{p_{k}_{i}})E(ε_{p_{k}_{j}}) = 0 , p_{k}_{i} ≢ p_{k}_{j} ,
E(ε_{p_{k}_{i}}ε_{p_{k}_{j}}) = σ_{p_{k}_{i}}^{2} , p_{k}_{i} ≡ p_{k}_{j} ,

∀ i, j, k : i ≠ j ∧ x_{i} ≢ x_{j} ∧ p_{k}_{i} ≢ p_{k}_{j}

v_{ij} = 0 ,

∀ i, j, k : i ≠ j ∧ x_{i} ≡ x_{j} ∧ p_{k}_{i} ≡ p_{k}_{j}

v_{ij} = (\frac{∂y}{∂x})_{x=x_{i},y=y_{i}} (\frac{∂y}{∂x})_{x=x_{j},y=y_{j}} σ_{x_{i}}^{2} + ∑_{k} (\frac{∂y}{∂p_{k}} )_{y=y_{i},p_{k}=p_{k}_{i}} (\frac{∂y}{∂p_{k}} )_{y=y_{j},p_{k}=p_{k}_{j}} σ_{p_{k}_{i}}^{2} .

If the “independent” measuring error ε_{y} is a predominant part of the true regression error ε, the distribution of the random vector ε is close to normal and the maximum likelihood estimate of the unknown parameters is close to the solution of the over-determined set of m equations (for n unknown variables) minimizing the generalized least squares method criterion§.

All the elements of the vector ϑ of all the calibration functions implemented so far are “linear”, and thus the OPIUM program \underline{does not require} their initial estimates. When selected elements of the vector ϑ are the only unknown parameters, then a calculation is made of their unbiased estimates by the iterative

¶ indexing of variables was taken into account, i.e.

ε_{i} = (\frac{∂y}{∂x})_{x=x_{i},y=y_{i}} ε_{x_{i}} + ε_{y_{i}} + ∑_{k} (\frac{∂y}{∂p_{k}} )_{y=y_{i},p_{k}=p_{k}_{i}} ε_{p_{k}_{i}}

§ In addition to a point estimate of mean of estimate of the unknown parameter(s) the OPIUM program also offers classical analysis of this estimate in a linearization approximation (variance, t- and F -interval estimate of mean, significance test, correlation matrix, partial and multiple correlation coefficients etc.) as well as jack-knife and likelihood ratio analysis (interval estimate of mean, significance test). At each point, many types of residuals are calculated along with a number of scalar and vector quantities for indication of golden points, outliers and high leverage points. The OPIUM program also tests several hypotheses on homogeneity, independence and normality of the standardized residuals.

linear least squares method (the variance-covariance matrix V is in general a function of the unknown parameters), where biased estimates are also available, obtained by the ridge regression method.

1.4 Principal component analysis

The dependence of the measured variable y on the free concentrations c of components in a common spectrophotometric or “fast exchange” NMR experiment¶ can be expressed in a matrix form

A = ε C , A = δΠC/C †,
A ∼ (m, n) , ε ∼ (m, k) , δ ∼ (m, k) , C ∼ (k, n) ,

m is the number of wavelengths (NMR active nuclei), n is the number of solutions, k is the number of independent (“active”) components, A is the absorbance matrix (the matrix of chemical shifts) (a_{ij} ≡ y_{(i−1)×n+j}), ε is the matrix of molar absorption coefficients, δ is the matrix of the chemical shifts of pure components and C is the matrix of the free concentrations of components (a column per solution).

Obviously r(A) ≤ min(m, n, k). Thus, if both m > k and n > k, the true rank of the absorbance matrix (the matrix of chemical shifts) gives us a reasonable estimate of the number of light absorbing (δ influencing) components in a solution k ≥ r(A). Moreover, the number of observations can dramatically be reduced from m × n to k × n by means of the singular value decomposition of the matrix A (no more than k suitable linear combinations of the m original wavelengths are needed).

For m ≥ n

A = U_{n}^{T}ΣV ,
U_{n} ∼ (n, m) , Σ ∼ (n, n) , V ∼ (n, n) ,
U_{n}U_{n}^{T} = VV^{T} = 1_{n} ,
Σ = diag(σ_{1}, …, σ_{n}) , σ_{1} ≥ … ≥ σ_{n} ≥ 0 ,
ΣV = U_{n}εC , ΣV = U_{n}δΠC/C .

Since r(A) ≐ k, Σ ≐ diag(σ_{1}, …, σ_{k}, 0_{k+1}, …, 0_{n})§. Then

Σ_{k}V_{k} ≐ U_{k}εC , Σ_{k}V_{k} ≐ U_{k}δΠC/C ,

¶ see the definitions of calibration functions of types 7 and −5 † Since only the “active” components are taken into account, matrix Π is diagonal containig only the \underline{non-zero} diagonal elements of the matrix P (see the definition of a calibration function of type −5), i.e.

Π ∼ (k, k), r(Π) = r(P) = k.

§ To estimate the “effective” rank k of the matrix A (both the equations would hold exactly only if

where the matrices Σ_{k} ∼ (k, k), V_{k} ∼ (k, n) and U_{k} ∼ (k, m) are appropriate submatrices of the matrices Σ, V and U_{n}¶.

For m < n

A = U^{T}ΣV_{m} ,
U ∼ (m, m) , Σ ∼ (m, m) , V_{m} ∼ (m, n) ,
UU^{T} = V_{m}V_{m}^{T} = 1_{m} ,
Σ = diag(σ_{1}, …, σ_{m}) , σ_{1} ≥ … ≥ σ_{m} ≥ 0 ,
ΣV_{m} = UεC , ΣV_{m} = UδΠC/C .

Since r(A) ≐ k, Σ ≐ diag(σ_{1}, …, σ_{k}, 0_{k+1}, …, 0_{m}). Then

Σ_{k}V_{k} ≐ U_{k}εC , Σ_{k}V_{k} ≐ U_{k}δΠC/C ,

where the matrices Σ_{k} ∼ (k, k), V_{k} ∼ (k, n) and U_{k} ∼ (k, m) are appropriate submatrices of the matrices Σ, V_{m} and U.

1.5 �Trial and error� method

Let p, k, l are \underline{natural} numbers and m, n are \underline{non-negative integers} (m+n < p ∧ k, l > m) with the following meaning:

p  …  number of “unknown” complexes,
m  …  number of “unknown” complexes with higher priority (\underline{all} of

them must definitely be in every generated combination),

n … number of “unknown” complexes with lower priority, p − m − n … number of “unknown” complexes with common priority (at

least one of them must be in each generated combination),

k, l … lower (upper) limit of order of combination.

Then the number of combinations (chemical model hypotheses) is

\begin{matrix}n < k − m : & ∑{l−m \atop i=k−m} \binom{p−m}{i} \\ k − m ≤ n ≤ l − m : & ∑_{i=k−m}^{l−m} \binom{p−m}{i} − ∑_{i=k−m}^{n} \binom{n}{i} \\ n > l − m : & ∑_{i=k−m}^{l−m} \binom{p−m}{i} − ∑_{i=k−m}^{l−m} \binom{n}{i} \end{matrix} .

s_{i} = \sqrt{\frac{tr (Σ^{2}) − ∑_{j=1}^{i} σ_{j}^{2}}{p − i}} , i = 0, …, p − 1 , p = min(m, n) .

there were no measuring errors), sequence of residual standard deviations is calculated^{3} Then inequalities s_{0} > … > s_{k−1} > σ_{y} ≥ s_{k} ≥ … ≥ s_{p} ≥ 0 should be valid, where σ_{y} is the standard deviation of the random variable y (absorbance, chemical shift).

¶ i.e. A’ ≐ ε’C, A’ ≐ δ’ΠC/C, where the matrices A’ = Σ_{k}V_{k}^{T} ∼ (k, n), ε’ = U_{k}^{T}ε ∼ (k, k) and δ’ = U_{k}^{T}δ ∼ (k, k) are “reduced” equivalents of the matrices A, ε and δ

For each combination, the OPIUM program creates a shortened information about the goodness of fit (containing quantities such as sum of squares of residuals (U), Akaike information criterion (AIC), Schwarz^{4}-Rissanen^{5}, Hannan- -Quinn^{6} and Anděl-Perez-Negrao^{7} statistic) and stores it into a file¶.

2 Data Analysis †§

2.1 Potentiometry

2-methylaminomethylpyrrole-4-sulphonic acid (H_{2}psma), Ag^{+} The 1st example^{8} is the potentiometric titration of 2-methylaminomethyl-

pyrrole-4-sulphonic acid and Ag^{+} (of known analytical concentration c_{1} of H^{+}, c_{2} of psma2− and c_{3} of Ag^{+}, the last two at the variable molar ratio 1 : 1 up to 2 : 1) with potassium hydroxide (of known analytical concentration c’_{1} < 0 of H^{+}). The glass electrode calibrated in concentration scale (pH = − log c_{1}) was used and the dependence of − log c_{1} on the added volume of the titrant was measured‡. The origin of individual parameters is displayed as an “extension” of the matrix q (the chemical model is taken from Stünzi H. et al^{8}).

H^{+} psma2− Ag^{+}
1  0  0  H^{+}  c_{1}, c’_{1}
0  1  0  psma2−  c_{2}
0  0  1  Ag^{+}  c_{3}
−1  0  0  OH−  β_{−100}
−1  0  1  AgOH  β_{−101}
1  1  0  Hpsma−  β_{110}
2  1  0  H_{2}psma  β_{210}
1  1  1  AgHpsma  β_{111}
2  2  1  AgH_{2}psma_{2}−  β_{221}
0  2  2  Ag_{2}psma_{2} 2−  β_{022}
1  3  3  Ag_{3}Hpsma_{3} 2−  β_{133}
0  3  3  Ag_{3}psma_{3} 3−  β_{033}
quantity  unit
mass amount  mmol
volume  ml
concentration mol · l^{−1}

¶ quite suitable for processing with standard Unix tools (grep, sort, awk etc.)

† common units are used in the rest of the text (and will no more be mentioned)

§ Unit variance-covariance matrix V = 1 is used in all the calculations to enable the comparison of the estimates found by different programs. On this assumption, all the programs using the same least squares criterion should give \underline{absolutely} the same results, i.e. point estimate(s) of mean(s) of the estimate(s) of the unknown parameter(s).

‡ 1 calibration function of type −1 is used

The results are presented in Table 1.

HYPERQUAD^{9} OPIUM
log β_{111}  18.00(1)  17.98(1)
log β_{221}  36.45(1)  36.43(1)
log β_{022}  20.17(1)  20.12(2)
log β_{133}  41.01(2)  41.00(2)
log β_{033}  31.39(3)  31.40(2)
10^{4} · U  115.4  115.4

Table 1.

When applying “trial and error” method, the ability of the program to get the estimates of the unknown parameters even from their bad initial estimates can save much time¶. That is why the robustness of both the programs used is tested before the best model is found. The results are presented in Table 2.

initial estimates  HYPERQUAD  OPIUM
20, 35, 20, 40, 30  +  +
10, 30, 10, 30, 20  no convergence^{a}  +
20, 35, 20, 0, 0  fails in 23rd iteration^{b}  +
0, 0, 0, 40, 30  +  +
0, 0, 0, 0, 0  no convergence^{a}  +
20, 35, 20, 50, 40  +  +
30, 30, 30, 40, 30 fails in 174th iteration^{b}  +
25, 30, 25, 35, 35  fails in 10th iteration^{b}  +
15, 40, 15, 45, 25  +  +
10, 45, 10, 50, 20  +  +
20, 20, 20, 20, 20  fails in 24th iteration^{b}  +
50, 50, 50, 50, 50  fails in 2nd iteration^{b}  +^{d}
0, 50, 0, 50, 0  fails in 6th iteration^{b}  +
50, 0, 50, 0, 50  fails in 1st iteration^{c}  +^{d}
50, 50, 50, 0, 0  fails in 1st iteration^{c}  +^{d}
0, 0, 0, 50, 50  +  +

a possibly a very very slow convergence b MS DOS◯^{c} collapses, rebooting is necessary c floating point error, program exits normally d takes approximately 5 minutes on i486DX 33MHz

Table 2.

¶ not to mention the possibility of passing over an existing solution

To determine the best model, we chose 23 “possible” complexes (corresponding rows of the matrix of stoichiometric coefficients q (H^{+}, psma2−, Ag^{+}) are

3. 2 2 1 6. 2 2 2  9. 2 3 2 12. 2 3 3 15. 2 4 2 18. 2 4 3 21. 2 4 4
1. 1 1 1 4. 1 2 1 7. 1 2 2 10. 1 3 2 13. 1 3 3 16. 1 4 2 19. 1 4 3 22. 1 4 4
2. 0 1 1 5. 0 2 1 8. 0 2 2 11. 0 3 2 14. 0 3 3 17. 0 4 2 20. 0 4 3 23. 0 4 4

and obtained ∑_{i=1}^{5} \binom{23}{i} = 44551 combinations of order 1 up to 5¶. To demonstrate the power of the OPIUM program, a “brute force” approach

was used and \underline{all the combinations} were tested†. 20818 solutions were found, but only 4035 of them are “positive” (i.e. point estimates of all the unknown stability constants are positive)§.

1. {5} o.k. 1.1539e-2 -1568.7135 -9.3825 -9.4384 1.7273e-4 "agpsma +[1 3 8 13 14]"
2. {5} o.k. 1.6055e-2 -1514.2158 -9.0522 -9.1081 2.4033e-4 "agpsma +[1 3 8 14 22]"
3. {5} o.k. 1.6485e-2 -1509.8532 -9.0257 -9.0816 2.4677e-4 "agpsma +[1 3 8 13 23]"
4. {5} o.k. 1.8080e-2 -1494.6209 -8.9334 -8.9893 2.7064e-4 "agpsma +[1 3 8 22 23]"
5. {5} o.k. 2.3136e-2 -1453.9293 -8.6868 -8.7427 3.4633e-4 "agpsma +[1 3 8 14 18]"

. .

13. {4} o.k. 3.7961e-2 -1374.2290 -8.2288 -8.2736 4.9893e-4 "agpsma +[1 3 8 22]"

. .

19. {4} o.k. 4.4369e-2 -1348.4928 -8.0728 -8.1176 5.8315e-4 "agpsma +[1 3 8 13]"

. .

4035. {1} o.k. 150.8192 -12.8275 -5.2839e-2 -6.4022e-2 1.1762 "agpsma +[6]"

As you can see, the best model is really the one analyzed above. For comparison, several “worse” alternatives are also presented in Table 3.

Some attempts to improve the goodness of fit and to confirm the determined model (q, β) by the ESI approach^{2} are presented in Table 4. The linear parameter E_{0}‡ represents the systematic error of the pH measurements here and should be 0, the original values of c_{2} (the analytical concentrations of psma2−) are 0.01001, 0.0203, 0.005, 0.00202, 0.01005, 0.00303, 0.00404¶, the correct values of q_{53}, q_{63}, q_{73} (the stoichiometric coefficients of Ag^{+} in the first three “unknown” complexes, i.e. 0 1 1, 2 2 1, 0 2 2) are 1, 1, 2.

¶ if some preliminary analysis of the data had been done, the number of combinations could have been reduced substantially^{8} † the task was run overnight on Alpha 21164A 300MHz, which is about 20 × faster than i486DX 33MHz § each line in the table gives a shortened information about respective model and goodness of fit: 5 is the order of the combination (i.e. the number of unknown complexes), o.k. means a positive solution ('N' would mean a negative solution), sum of squares of residuals (U), Akaike information criterion (AIC), Schwarz-Rissanen, Hannan-Quinn (1), Anděl-Perez-Negrao (1, \frac{1}{4}) statistic, agpsma is the name of the input file being processed and +[1 3 8 13 14] is the set of complexes forming the model (the quoted information is in the format required by the OPIUM program as the command line parameter to “rerun” the model and produce the complete output file) ‡ 7 calibration functions of type 2 are used, S = −1 ¶ the E_{0} and c_{2} parameters are supposed to be individual properties of each of the 7 titration curves

1.  2.  3.  4.  5.  13.  19.
log β_{111}  17.98(1) 17.93(1) 18.00(1)  17.95(1)  18.00(1)  17.96(2)  18.00(2)
log β_{221}  36.43(1) 36.50(1) 36.40(1)  36.48(1)  36.43(1)  36.46(1)  36.39(2)
log β_{022}  20.12(2) 20.17(2) 20.23(1) 20.246(9) 20.13(2) 20.339(6) 20.350(7)
log β_{133}  41.00(2)  -  41.02(2)  -  -  -  40.78(3)
log β_{033}  31.40(2) 31.27(3)  -  -  31.39(4)  -  -
log β_{243}  -  -  -  -  59.37(3)  -  -
log β_{144}  -  52.39(2)  -  52.42(2)  -  52.21(3)  -
log β_{044}  -  -  42.55(3)  42.41(4)  -  -  -
10^{4} · U  115.4  160.6  164.9  180.8  231.4  379.6  443.7
AIC  −1569  −1514  −1510  −1495  −1454  −1374  −1348

Table 3.

E_{0}, β  c, β  q, β  E_{0}, c, β  c, q, β  boundary condition
E_{0}  −0.006(6)  ·  ·  0.005(5)  ·  -
E_{0}  0.008(7)  ·  ·  −0.005(5)  ·  -
E_{0}  0.007(6)  ·  ·  0.002(5)  ·  -
E_{0}  0.007(5)  ·  ·  −0.007(4)  ·  -
E_{0}  0.007(4)  ·  ·  0.012(5)  ·  -
E_{0}  0.015(5)  ·  ·  0.010(4)  ·  -
E_{0}  0.010(4)  ·  ·  0.019(4)  ·  -
10^{2} · c_{2}  ·  1.006(1)  ·  1.005(1)  1.003(2)  c_{1} = 0.01967 + 2c_{2}
10^{2} · c_{2}  ·  2.003(3)  ·  1.991(2)  1.998(4)  c_{1} = 0.01987 + 2c_{2}
10^{3} · c_{2}  ·  4.974(7)  ·  4.963(7)  4.96(1)  c_{1} = 0.01975 + 2c_{2}
10^{3} · c_{2}  ·  2.018(4)  ·  1.991(5)  2.014(7)  c_{1} = 0.00398 + 2c_{2}
10^{2} · c_{2}  ·  0.999(1)  ·  1.007(2)  0.998(3)  c_{1} = 0.01970 + 2c_{2}
10^{3} · c_{2}  ·  2.996(5)  ·  3.015(7)  2.99(1)  c_{1} = 0.00398 + 2c_{2}
10^{3} · c_{2}  ·  4.015(6)  ·  4.09(1)  4.01(2)  c_{1} = 0.00398 + 2c_{2}
q_{53}  ·  ·  0.99(2)  ·  0.99(1)  -
q_{63}  ·  ·  1.06(2)  ·  1.039(6)  -
q_{73}  ·  ·  2.048(9)  ·  2.01(2)  -
log β_{111}  17.99(1)  17.984(4) 17.94(4)  17.963(8)  17.95(2)  -
log β_{221}  36.44(1)  36.442(5) 36.58(4)  36.48(1)  36.55(2)  -
log β_{022}  20.21(3)  20.11(2)  20.31(4)  20.13(1)  20.13(5)  -
log β_{133}  40.97(4)  40.996(9) 40.95(2)  40.97(3)  41.00(1)  -
log β_{033}  31.30(7)  31.49(3)  31.36(3)  31.51(5)  31.53(2)  -
10^{4} · U  50.80  17.84  94.76  6.018  12.29
AIC  −1690  −1863  −1595  −2028  −1918

· means the original (correct) value

Table 4.

2.2 Spectrophotometry

methionine (Hmet), SCN−, Pd^{2+}

The 2nd example^{10,11} is the spectrophotometric “titration” of methionine¶, SCN− and Pd^{2+} at the fixed molar ratio 2 : 1000 : 1 and pH 0 up to 11. The free concentration c_{1} of H^{+} (i.e. pH) and the analytical concentration C_{2} of met−, C_{3} of SCN− and C_{4} of Pd^{2+} are known in each of the 30 solutions. The dependence of the absorbance at 25 wavelengths (260 up to 380 nm) on the free concentration c_{1} (for fixed analytical concentrations C_{2}, C_{3} and C_{4}) was measured†. The origin of individual parameters is displayed as an “extension” of the matrix q (the chemical model is taken from Stypinski B.^{10}).

H^{+} met− SCN− Pd^{2+}
1  0  0  0  H^{+}  30 × c_{1}
0  1  0  0  met−  30 × C_{2}
0  0  1  0  SCN−  30 × C_{3}
0  0  0  1  Pd^{2+}  30 × C_{4}
1  1  0  0  Hmet  β_{1100}
2  1  0  0  H_{2}met^{+}  β_{2100}
0  0  4  1  Pd(SCN)_{4} 2−  25 × a_{3}, β_{0041}
0  1  2  1  Pdmet(SCN)_{2}−  25 × a_{4}, β_{0121}
0  2  0  1  Pdmet_{2}  25 × a_{5}, β_{0201}

Principal component analysis gives the set of residual standard deviations of absorbance 0.2581, 0.0466, 0.035, 0.0009, i.e. k ≈ 3 (there are at least 3 light absorbing species). The results are presented in Table 5.

SPECFIT^{12}  OPIUM
log β_{0121}  31.561(7)  31.561(8)
log β_{0201}  34.33(2)  34.33(3)
10^{4} · U  98.96  104.9

Table 5.

When only 3 principal components of the absorbance matrix are selected, the OPIUM program gives the sum of squares of residuals 10^{4} · U = 98.96 (which results in smaller values of the estimated standard deviations).

2-methyl-3-hydroxy-4H-pyran-4-one (maltol, maltol), UO_{2}^{2+} The 3rd example^{13} is the spectrophotometric “titration” of maltol and UO_{2}^{2+} at the variable molar ratio 10 : 1 up to 300 : 1 and pH 1 up to 9. The free concentration c_{1} of H^{+} (i.e. pH) and the analytical concentration C_{2} of maltol and C_{3} of UO_{2}^{2+} are known in each of the 148 solutions. The dependence of the absorbance at 3 wavelengths (390, 400 and 410 nm using 2 different optical lengths) on the free concentration c_{1} and the analytical concentrations C_{2} and C_{3} was measured¶. The origin of individual parameters is displayed as an “extension” of the matrix q (the chemical model is taken from Havel J.^{13}).

¶ 2-amino-4-(methylthio)butanoic acid † 25 calibration functions of type 7 are used

H^{+} maltol UO_{2}^{2+}
1  0  0  H^{+}  148 × c_{1}
0  1  0  maltol  148 × C_{2}
0  0  1  UO_{2}^{2+}  148 × C_{3}
−1  1  0  H−1maltol−  β_{−110}
0  1  1  UO_{2}maltol^{2+}  3 × a_{2}, β_{011}
−1  1  1  UO_{2}H−1maltol^{+}  3 × a_{3}, β_{−111}
−1  2  1  UO_{2}H−1maltol_{2}^{+}  3 × a_{4}, β_{−111}
−2  2  1  UO_{2}H−2maltol_{2}  3 × a_{5}, β_{−221}
−2  3  1  UO_{2}H−2maltol_{3}  3 × a_{6}, β_{−231}
−3  3  1  UO_{2}H−3maltol_{3}−  3 × a_{7}, β_{−331}
−3  1  1  UO_{2}H−3maltol−  3 × a_{8}, β_{−311}
−6  4  2  (UO_{2})_{2}H−6maltol_{4} 2−  3 × a_{9}, β_{−642}
−6  3  2  (UO_{2})_{2}H−6maltol_{3} 2−  3 × a_{10}, β_{−632}
−6  1  2  (UO_{2})_{2}H−6maltol2−  3 × a_{11}, β_{−612}

In spite of the enormous number of solutions, the absorbance data do not contain enough information for computation of estimates of all the unknown parameters (a, β) together (and no “overall” calculation has ever been performed). The 10 unknown complexes were determined by dividing the studied pH range into several parts and treating them separately.

We did not follow this procedure, but searched for the best “overall” model. For the sake of simplicity, we chose the 10 already “determined” complexes and obtained ∑_{i=1}^{10} \binom{10}{i} = 1023 combinations of order 1 up to 10.

There was no problem to test \underline{all the combinations}†. 740 solutions were found, 535 of them are “positive” (i.e. point estimates of all the unknown stability constants are positive), but most of them are highly correlated. That is why only positive solutions of orders 1 up to 5 are listed (421 solutions).

1. {5} o.k. 1.5468e-2 -4517.5667 -9.9441 -10.0559 1.9543e-4 "uo2malt +10 9 5 3 1"
2. {5} o.k. 1.6792e-2 -4481.1064 -9.8620 -9.9737 2.1216e-4 "uo2malt +10 7 5 3 1"
3. {5} o.k. 1.7515e-2 -4462.3811 -9.8198 -9.9316 2.2130e-4 "uo2malt +8 7 5 3 1"

¶ 3 calibration functions of type 8 are used † the task took several minutes on R10000 175MHz, which is about 20 × faster than i486DX 33MHz

4. {5} o.k. 1.7573e-2 -4460.9123 -9.8165 -9.9283 2.2203e-4 "uo2malt +10 8 5 3 1"
5. {5} o.k. 2.0239e-2 -4398.2185 -9.6753 -9.7870 2.5570e-4 "uo2malt +9 7 5 3 1"
6. {5} o.k. 2.1360e-2 -4374.2718 -9.6214 -9.7331 2.6987e-4 "uo2malt +9 8 5 3 1"
7. {5} o.k. 2.6417e-2 -4279.9353 -9.4089 -9.5206 3.3376e-4 "uo2malt +10 9 7 5 3"
8. {4} o.k. 2.7340e-2 -4272.6845 -9.4389 -9.5283 2.8653e-4 "uo2malt +[2 4 6 8]"
9. {5} o.k. 3.3617e-2 -4172.9169 -9.1679 -9.2796 4.2473e-4 "uo2malt +10 9 8 3 1"
10. {5} o.k. 3.4558e-2 -4160.6572 -9.1403 -9.2520 4.3662e-4 "uo2malt +10 9 8 5 3"
11. {4} o.k. 3.5413e-2 -4157.8077 -9.1801 -9.2695 3.7114e-4 "uo2malt +[2 4 6 7]"

. .

421. {1} o.k. 12.6952 -1570.2417 -3.4906 -3.5130 5.3995e-2 "uo2malt +[1]"

The best models (of orders 5 and 4) are presented in Table 6.

1.  2.  3.  5.  8.  11.
log β_{−111}  −0.40(2)  −0.40(2)  −0.39(2)  −0.39(2)  −0.38(2)  −0.40(3)
log β_{−221}  −2.33(3)  −2.31(4)  −2.30(4)  −2.28(4)  −2.25(5)  −2.33(5)
log β_{−331}  −7.57(4)  −7.59(4)  −7.49(4)  −7.50(5)  −7.42(5)  −7.48(6)
log β_{−311}  −12.20(6)  -  -  -  -  −11.92(6)
log β_{−642}  −15.49(8) −15.1(1)  -  −15.45(9) −15.6(1)  -
log β_{−632}  -  −17.4(1) −17.17(9)  -  -  -
log β_{−612}  -  -  −22.7(1)  −23.8(1)  -  -
10^{4} · U  154.7  167.9  175.2  202.4  273.4  354.1
AIC  −4518  −4481  −4462  −4398  −4273  −4158

Table 6.

2.3 �Fast exchange� NMR

4,9-diaza-1,12-dodecanediamine (spermine, sprm)

The 4th example^{14,15} is the ^{1}H+^{13}C NMR “titration” of spermine. The free concentration c_{1} of H^{+} (i.e. pH) and the analytical concentration C_{2} of sprm are known in each of the 20 solutions. There are 10 different types of nuclei (5×H, 5×C) in the molecule and thus the dependence of 10 chemical shifts on the free concentration c_{1} was measured¶. The origin of individual parameters is displayed as an “extension” of the matrix q.

H^{+} sprm
1  0  H^{+}  20 × c_{1}
0  1  sprm  10 × a_{−2}, 20 × C_{2}
1  1  Hsprm^{+}  10 × a_{1}, β_{11}
2  1  H_{2}sprm^{2+}  10 × a_{2}, β_{21}
3  1  H_{3}sprm^{3+}  10 × a_{3}, β_{31}
4  1  H_{4}sprm^{4+}  10 × a_{4}, β_{41}

Principal component analysis gives the set of residual standard deviations

¶ 10 calibration functions of type −5 are used

of chemical shift 28.51, 0.8013, 0.0522, 0.0252, 0.0017, 0.0009, i.e. k ≈ 5 (there are at least 5 δ influencing species, but even for k = 1 reliable estimates can be obtained). The results are presented in Table 8.

HYPNMR^{14,15} OPIUM
log β_{11}  11.20(2)  11.20(2)
log β_{21}  21.49(6)  21.49(6)
log β_{31}  30.53(8)  30.53(8)
log β_{41}  38.7(1)  38.7(1)
10^{4} · U  37.09  37.09

Table 8.

The NMR chemical shift(s) of a pure weak acid in a buffered solution is independent on the concentration of the acid¶. That is why no of the 20 analytical concentrations C_{2} can be “marked” unknown. Formally, it is possible to calculate any of the 20 free concentrations c_{1} of H^{+}, but it can hardly be of any practical importance. Moreover, only a few of the free concentrations can be calculated because of the strong correlation of all the estimates†.

2-aminoethylphosphonic acid (aep2−)

The 5th example^{16,17} is the ^{31}P NMR titration of 2-aminoethylphosphonic acid (of known mass amount n_{1} of H^{+} and n_{2} of aep2−) with tetramethylammonium hydroxide (of known analytical concentration c’_{1} < 0 of H^{+}). There is one phosphorus nucleus in the molecule and thus the dependence of one chemical shift on the added volume of the titrant was measured§. The origin of individual parameters is displayed as an “extension” of the matrix q.

H^{+} aep2−
1  0  H^{+}  n_{1}, c’_{1}
0  1  aep2−  a_{−2}, n_{2}
−1  0  OH−  β_{−10}
1  1  Haep−  a_{2}, β_{11}
2  1  H_{2}aep  a_{3}, β_{21}
3  1  H_{3}aep^{+}  a_{4}, β_{31}

The hypersurface of the sum of squares of residuals with respect to the parameters a, β is “flat” (in the 8-dimensional space) because of the strong correlation of the estimates of parameters a_{4}, β_{21} and β_{31}. At least one of the parameters n, V, c’ \underline{has to be} “marked” unknown to make the problem solvable. Several choices¶ are presented in Table 9. The linear parameters a_{−2}, a_{2}, a_{3} and a_{4} are the^{31}P NMR chemical shifts of pure components aep2−, Haep−, H_{2}aep and H_{3}aep^{+}, the original values of n_{1} and n_{2} (the mass amounts of H^{+} and aep2− in the titrand) are 2.6371 and 0.8760, the original value of c’_{1} (the analytical concentration of H^{+} in the titrant) is −0.09912.

¶ all the analytical concentrations C_{2} ≡ C in the input file could be replaced by random numbers and the results, i.e. estimates of the unknown parameters, would be the same † the number of solutions must be substantially greater than the number of unknown parameters § 1 calibration function of type −5 is used

a, n_{1}, β  a, n_{2}, β  a, c’_{1}, β  a, n, β  a, n_{1}, c’_{1}, β  a, n_{2}, c’_{1}, β
a_{−2}  20.16(1)  20.18(1)  20.15(1)  20.17(1)  20.17(1)  20.17(1)
a_{2}  16.813(1) 16.812(1)  16.815(1)  16.812(1)  16.812(1)  16.812(1)
a_{3}  19.272(1) 19.263(2)  19.278(1)  19.265(2)  19.265(2)  19.265(2)
a_{4}  25(1)  23.4(5)  28(2)  23.7(6)  23.7(6)  23.7(6)
n_{1}  2.6430(3)  ·  ·  2.639(1)  2.651(2)  ·
n_{2}  ·  0.8708(3)  ·  0.872(1)  ·  0.8715(7)
c’_{1}  ·  ·  −0.09879(2)  ·  −0.0996(1)  −0.09906(5)
log β_{11}  11.361(6) 11.370(5)  11.356(6)  11.367(6)  11.369(6)  11.367(6)
log β_{21}  18.71(5)  18.77(5)  18.66(6)  18.76(5)  18.76(5)  18.75(5)
log β_{31}  19.7(1)  19.95(8)  19.4(1)  19.9(1)  19.9(1)  19.9(1)
10^{4} · U  16.31  14.92  19.68  14.76  14.77  14.76
AIC  −1414  −1426  −1391  −1425  −1425  −1425

· means the original value

Table 9.

2.4 �Slow exchange� NMR

1-(2-aminoethylamino)-1-deoxy-D-galacticol (galac), Cd^{2+} The 6th example^{18} is the ^{113}Cd NMR “titration” of 1-(2-aminoethylamino)-

-1-deoxy-D-galacticol and Cd^{2+} at the variable molar ratio 1 : 5 up to 2 : 1 and approximately constant pH ≈ 6.5. The free concentration c_{1} of H^{+} (i.e. pH) and the analytical concentration C_{2} of galac and C_{3} of Cd^{2+} are known in each of the 9 solutions. There are 3 peaks in the NMR spectrum (assigned to Cd^{2+}, Cdgalac^{2+} and Cdgalac_{2}^{2+}) and thus the dependence of three \underline{normalized} integral intensities on the free concentration c_{1} and the analytical concentrations C_{2} and C_{3} was measured†. The origin of individual parameters is displayed as an “extension” of the matrix q.

H^{+} galac Cd^{2+}
1  0  0  H^{+}  9 × c_{1}

¶ some of them are a bit obscure, probably the best approach should be the calculation of n_{2} using the “natural” boundary condition for n_{1}, but the information was not available † 1 calibration function of type −6 is used (with 3 different specifiers i, see the definition of the function)

H^{+} galac Cd^{2+}
0  1  0  galac  9 × C_{2}
0  0  1  Cd^{2+}  9 × C_{2}
1  1  0  Hgalac^{+}  β_{110}
2  1  0  H_{2}galac^{2+}  β_{210}
0  1  1  Cdgalac^{2+}  β_{011}
0  2  1  Cdgalac_{2}^{2+}  β_{021}

The results are presented in Table 10.

Huskens J. et al^{18} OPIUM
log β_{011}  4.6(3)  4.63(6)
log β_{021}  8.7(4)  8.6(1)
10^{4} · U  658.4

Table 10.

113Cd NMR integral intensity data are usually of a poor quality. On the other hand, the chemical shifts of individual peaks give us reliable information about both the number of “active” components and their structure. No matter how imprecise the integral intensities are, they should agree with the proposed chemical model. Such an agreement is tested by the ESI approach and the results are presented in Table 11.¶† The correct values of q_{32}, q_{33}, q_{42}, q_{43} (the stoichiometric coefficients of galac and Cd^{2+} in both the unknown complexes, i.e. 0 1 1, 0 2 1) are 1, 1, 2, 1.

q_{32}, β  q_{33}, β  q_{42}, β  q_{43}, β  q_{32}, q_{42}, β  q_{33}, q_{43}, β
q_{32}  0.83(7)  ·  ·  ·  0.88(8)  ·
q_{33}  ·  1.3(1)  ·  ·  ·  1.3(1)
q_{42}  ·  ·  3.8(6)  ·  2.8(5)  ·
q_{43}  ·  ·  ·  0.7(2)  ·  0.7(2)
log β_{011}  3.8(4)  4.8(1) 5.0(2) 4.8(2)  4.1(5)  4.9(2)
log β_{021}  8.4(1)  8.4(1)  18(3)  8.7(2)  13(3)  8.4(2)
10^{4} · U  528.4  535.4  496.9  634.3  485.1  502.7

· means the correct value

Table 11.

As you can see, the correlation of the estimates of the parameters β and q makes the analysis almost impossible.

¶ 1 calibration function of type −6 is used † if the calibration function of type −4 were used (α = 1, specifier d set to the reciprocal analytical concentration C_{3} of Cd^{2+} in each of the 9 solutions d = 1/C_{3} ≡ 1/C), the estimate of any of the stoichiometric coefficients of Cd^{2+} would be incorrect (C is a function of q)

\underline{3 References}

1. Kývala M., Lukeš I., CHEMOMETRICS ’95, Pardubice, Czech Republic, 1995, Book of Abstracts p. 63

2. Havel J., Meloun M., Talanta 33 (1986) 525

3. Kankare J. J., Anal. Chem. 42 (1970) 1322

4. Schwarz G., Ann. Statist. 6 (1978) 461

5. Rissanen J., Automatica 14 (1978) 465

3. Hannan E. J., Quinn B. G., J. Roy. Statist. Soc. B41 (1979) 190

7. Anděl J., Perez M. G., Negrao A. I., Kybernetika 17 (1981) 514

8. Stünzi H., Anderegg G., J. Coord. Chem. 7 (1978) 239

9. Gans P., Sabatini A., Vacca A., Talanta 43 (1996) 1739

10. Stypinski B., Dissertation ETH Nr. 11713, Zürich, 1996

11. Anderegg G., personal communication

12. Gampp H., Maeder M., Meyer C. J., Zuberbühler A. D., Talanta 32 (1985) 257

13. Havel J., Pure Appl. Chem. 34 (1972) 369

14. Frassineti C., Ghelli S., Gans P., Sabatini A., Moruzzi M. S., Vacca A., Anal. Biochem. 231 (1995) 374

15. http://www.chim1.unifi.it/group/vacsab/hypnmr.htm

16. Ollig J., Hägele G., Computers Chem. 19 (1995) 287

17. Hägele G., personal communication

18. Huskens J., Lammers H., Bekkum H. van, Peters J. A., Magn. Reson. Chem. 32 (1994) 691

19. http://www.natur.cuni.cz/~kyvala/opium.html
20. http://www.chemi.muni.cz/~lubal