Description of the OPIUM Computer Program

The main feature of the OPIUM program^{1} (besides extensive analysis of both the estimates and the residuals) is a strict separation of the dependence of a response variable y (electromotive force of a cell in potentiometry, absorbance in UV/VIS spectrophotometry, chemical shift in “fast exchange” NMR spectroscopy, integral intensity of a peak in “slow exchange” NMR spectroscopy, reaction heat in enthalpimetry, difference of resistances in VPO, effective mobility in capillary zone electrophoresis of a labile system, conductivity in conductimetry etc.) on the calibration and the 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 solution, V is the volume of the solution, c’_{j} is the analytical concentration of the j-th reactant in the burette and x is the added volume (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, any subset of the parameters ϑ, c, n, V, c’, q, β can be “marked” unknown and calculated. Moreover, simple boundary conditions (linear constraints) 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 estimates of the unknown parameters. 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).

The OPIUM program is robust enough to usually find a sole global minimum or at least one of local minima even from roughly approximate initial estimates of the unknown stability constants† (i.e. without any human assistance). This feature allows us to test 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 under any operating system). 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. If the chemical model is sufficiently simple and the data are sufficiently precise, it is possible to attempt a direct computation of all the unknown stability constants and all the unknown stoichiometric coefficients without prior testing of hundreds or even thousands of hypotheses.

¶ a_{1}, a_{2}, a_{3} could be e.g. molar absorption coefficients

† 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, there is frequently just one local minimum on the hypersurface of the sum of squares of residuals and the method used for optimization of the unknown parameters usually ensures rapid convergence to the only solution. If they are not, some other local minima often exist and better initial estimates of β might be required to find the “correct” one. Fortunately, the initial estimates do not need to be exceedingly precise. The program can mostly get along any lower initial estimates of the unknown stability constants (owing to the problems of underflow we do not recommend to set the initial estimates of β_{i} close to the minimum positive real number representable by the computer being used, i.e. estimates of log β_{i} lesser than about −200 on computers using standard 8-byte real numbers). Initial estimates higher than the optimum values by several tens of orders of magnitude raise more serious problems. However, the program can usually find the correct solution even from those.

§ while being somewhere else and/or doing something more exciting

In the whole program, the first derivatives are calculated analytically, while the second derivatives (required in calculation of the maximum parametric and internal curvature of the model and Box bias of the estimate and residuals) are calculated numerically.

The percentage points of the selected distributions are calculated analytically directly from their definitions (where necessary, using the value of the error function erf x = \frac{2}{\sqrt{π}} ∫_{0}^{x} e^{−p^{2}}dp). The corresponding critical values are calculated by the Newton method.

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 function f with respect to the logarithms of the unknown free concentrations of reactants (elements of the symmetric Jacobi matrix) 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 .

The OPIUM program combines several methods of calculation of the unknown free concentrations of reactants. Besides the Newton-Raphson method with the restricted† shift vector δ = −J^{−1}f (using the Cholesky decomposition of the Jacobi matrix) and the Marquardt nonlinear least squares method with the restricted shift vector δ = −(J^{T}J + λ1)^{−1}J^{T}f (using the decomposition of the Jacobi matrix to its eigenvalues and eigenvectors)§ a robust (but slow) non-derivative method has been adopted^{3,4} (and slightly modified) with the restricted shift

C_{k} ≥ 0 : δ_{k} = (p_{k} + ϱ) [ln (C_{k} − B_{k}) − ln A_{k}] ,
C_{k} ≤ 0 : δ_{k} = (p_{k} + ϱ) (ln |B_{k}| − ln |C_{k} − A_{k}|) ,

A_{k} = c_{k} + ∑_{{i=1 \atop q_{ik}>0}}^{m} q_{ik}c_{n+i} , B_{k} = ∑_{{i=1 \atop q_{ik}<0}}^{m} q_{ik}c_{n+i} , A_{k} + B_{k} = C_{k} , k ∈ ⟨1, r⟩ ,

p_{k} + ϱ ∈ \{3, 2, 1, \frac{1}{2}, \frac{1}{3}, \frac{1}{4}, \frac{1}{5}, \frac{1}{6}, …\} .

¶ 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

† the elements of the shift vector δ must lie in the interval ⟨−1, 1⟩

§ the λ modifier of step length and direction is initialized with the value 0.02 tr (J^{T}J)/r

In each iteration, \underline{only one} of the free concentrations is corrected (mostly the one with the highest absolute value of the shift). While ϱ changes adaptively in the course of the computation, p_{k} is a characteristic constant of the k-th reactant. We choose it as the reciprocal value of one half of the weighted average of one (representing the reactant) and m absolute values of the elements of the k-th column of the matrix of stoichiometric coefficients q with the weights

(1 + ω)^{log c_{k}} , (1 + ω)^{log |q_{ik}|+log β_{i}+q_{ik} log c_{k}} , i = 1, …, m ,

ω is a conveniently selected empirical constant ω ∼ 0.001 − 0.01. Using mathematical symbols,

\frac{2}{p_{k}} = \frac{(1 + ω)^{log c_{k}} + ∑_{i=1}^{m} |q_{ik}|(1 + ω)^{log |q_{ik}|+log β_{i}+q_{ik} log c_{k}}}{(1 + ω)^{log c_{k}} + ∑_{i=1}^{m}(1 + ω)^{log |q_{ik}|+log β_{i}+q_{ik} log c_{k}}}

(the free concentration of the k-th reactant is replaced with its initial estimate).

Initial estimates of the unknown free concentrations of reactants are derived (at the first point of each titration section) using the known/fixed analytical concentrations of the reactants and the known/fixed ionic product of water (if it is defined). The logarithms of the unknown free concentrations are calculated with an accuracy to about one significant digit by the non- -derivative method similar to the algorithms used in the programs COMICS^{3} and COMPLEX^{4}. To ensure a smooth convergence of the method of optimization of the unknown parameters, we need to calculate the logarithms of the free concentrations with an accuracy to approximately 13 or 14 significant digits (on computers using standard 8-byte real numbers). The Newton-Raphson method usually takes care of this providing (close to the solution) a quadratic convergence rate. If it fails, the non-derivative method is applied again. If it is possible, we increase the accuracy of this computation as long as the ratio of the biggest and the smallest eigenvalue of the Jacobi matrix is greater than 10^{6}. However, only rarely the logarithms of the unknown free concentrations are calculated with an accuracy to more than 2 significant digits. From these estimates the Marquardt method usually finds the solution. If it does not, everything is repeated (with greater demands on the accuracy of the computation by the non- -derivative method and, if need be, more than once). This algorithm seems to be successful under all reasonable circumstances and mostly also very fast.

Thanks to a deliberate selection of the natural number k and of the rational number p_{k} + ϱ the non-derivative method is virtually always convergent and especially in the initial phase of the computation also relatively fast. The Newton-Raphson method ensures a quadratic convergence rate close to the solution. The Marquardt method is of use when the Jacobi matrix is almost singular and cannot safely be inverted even near the solution of the mass-balance equations.

The described algorithm obviously best resembles the algorithms used in the programs CHEMEQUIL-2^{5} and ES4EC^{6,7}.

At the second point of each titration section, the values computed at the first point are taken for initial estimates of the unknown free concentrations of reactants. The restricted linear extrapolation is applied at all the other points.

Starting with the second iteration, the values computed in the previous iteration are taken for initial estimates of the unknown free concentrations of reactants. If some of the unknown stability constants become negative, the simple non-derivative method cannot be applied even if both the derivative methods fail. Then it is necessary to multiply the negative stability constants with −1 and perform a restart of the method of optimization of the unknown parameters.

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} .

‡ excluding two exceptions (see the description of calibration functions of types 1 and −5)

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 ,

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.

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}}

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}}

¶ it is clear that

† 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).

§ including some “non-trivial” derivatives

ϑ ≡ 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}}

ϑ ≡ 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. Fe3+/Fe2+, Sn4+/Sn2+, Ce4+/Ce3+ 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}}

“Subtype” 9ln

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

i ≠ j ∈ ⟨1, n + m⟩

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

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

¶

p_{i} ∈ N+ \{0\} , 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 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).

\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

\underline{Type −5}

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

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

\frac{∂y}{∂c_{j}} = −a_{i}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}

A_{uv} = p_{v}c_{i} (a_{i} − 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 −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.

“Subtype” −5ln

y = − ln a_{i} − ln p_{i} − ln c_{i} + ln C ¶

\frac{∂y}{∂ ln a_{i}} = −1 , \frac{∂y}{∂c_{i}} = p_{i} (1 − \frac{C}{p_{i}c_{i}} ) /C = p_{i} (1 − e^{y+ln a_{i}}) /C

\frac{∂y}{∂c_{j}} = p_{j}/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}

A_{uv} = p_{v}c_{i} (1 − \frac{C}{p_{i}c_{i}} ) /C = p_{v}c_{i} (1 − e^{y+ln a_{i}}) /C , n + u = i , v = 1, …, n

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

Type −6

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

i ∈ ⟨1, n + m⟩, p_{i} ≠ 0

\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}

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

¶ p_{i} ≠ 0, ∀ j ≠ i : a_{j} = 0

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}).

“Subtype” −6ln

y = − ln p_{i} − ln c_{i} + ln C

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

\frac{∂y}{∂c_{j}} = p_{j}/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}

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

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

Type −7

y = α ∑_{i=1}^{n+m} c_{i}

ϑ ≡ α

\frac{∂y}{∂α} = ∑_{i=1}^{n+m} c_{i} , \frac{∂y}{∂c_{i}} = α , i = 1, …, n + m

Note: Function of type −7 describes the dependence of an \underline{extensive} linearly additive quantity y, usually a change of some colligative property (difference of resistances in VPO, osmotic pressure in membrane osmometry etc.) on the free concentrations c of all components.

Regression model

Let m is the number of points, i.e. the number of ordered pairs of measured quantities 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 the 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 ,
D(ε_{x}) = E(ε_{x}^{2}) = D(x) = σ_{x}^{2} ,
D(ε_{y}) = E(ε_{y}^{2}) = D(y) = σ_{y}^{2} ,
D(ε_{p_{k}}) = E(ε_{p_{k}}^{2}) = D(p_{k}) = σ_{p_{k}}^{2} ,
ε_{x} ∼ N(0, σ_{x}^{2}) , ε_{y} ∼ N(0, σ_{y}^{2}) , ε_{p_{k}} ∼ N(0, σ_{p_{k}}^{2}) .

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}

f(p, x) ≐ f(E(p), E(x)) + \frac{∂y}{∂x}ε_{x} + ∑_{k} \frac{∂y}{∂p_{k}} ε_{p_{k}}

(x = E(x) + ε_{x}, p_{k} = E(p_{k}) + ε_{p_{k}}) it follows that

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

E(ε) ≐ 0 ,

D(ε) ≐ E(ε^{2}) ≐ (\frac{∂y}{∂x})^{2} σ_{x}^{2} + σ_{y}^{2} + ∑_{k} (\frac{∂y}{∂p_{k}} )^{2} σ_{p_{k}}^{2} ,

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} ,

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

ε_{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}}

¶ all the necessary derivatives are given, except for

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

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} .

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

f(p, x) ≐ f(E(p), E(x)) + \frac{∂y}{∂x}ε_{x} + ∑_{k} \frac{∂y}{∂p_{k}} ε_{p_{k}} + \frac{1}{2} (\frac{∂^{2}y}{∂x^{2}} ε_{x}^{2} + ∑_{k} \frac{∂^{2}y}{∂p_{k}^{2}} ε_{p_{k}}^{2}) +

+ ∑_{k} \frac{∂^{2}y}{∂x∂p_{k}} ε_{x}ε_{p_{k}} + ∑_{k} ∑_{l>k} \frac{∂^{2}y}{∂p_{k}∂p_{l}} ε_{p_{k}}ε_{p_{l}}

(x = E(x) + ε_{x}, p_{k} = E(p_{k}) + ε_{p_{k}}) it follows that

ε ≐ \frac{∂y}{∂x}ε_{x} + ε_{y} + ∑_{k} \frac{∂y}{∂p_{k}} ε_{p_{k}} + \frac{1}{2} (\frac{∂^{2}y}{∂x^{2}} ε_{x}^{2} + ∑_{k} \frac{∂^{2}y}{∂p_{k}^{2}} ε_{p_{k}}^{2}) + ∑_{k} \frac{∂^{2}y}{∂x∂p_{k}} ε_{x}ε_{p_{k}} + ∑_{k} ∑_{l>k} \frac{∂^{2}y}{∂p_{k}∂p_{l}} ε_{p_{k}}ε_{p_{l}} ,

E(ε) ≐ \frac{1}{2} (\frac{∂^{2}y}{∂x^{2}} σ_{x}^{2} + ∑_{k} \frac{∂^{2}y}{∂p_{k}^{2}} σ_{p_{k}}^{2}) ,

E(ε) ≐ μ .

D(ε) = E(ε^{2}) − [E(ε)]^{2} ≐ (\frac{∂y}{∂x})^{2} σ_{x}^{2} + σ_{y}^{2} + ∑_{k} (\frac{∂y}{∂p_{k}} )^{2} σ_{p_{k}}^{2} + \frac{1}{4} \{(\frac{∂^{2}y}{∂x^{2}} )^{2} [E(ε_{x}^{4}) − σ_{x}^{4}] +

+ ∑_{k} (\frac{∂^{2}y}{∂p_{k}^{2}} )^{2} [E(ε_{p_{k}}^{4}) − σ_{p_{k}}^{4}]\} + ∑_{k} (\frac{∂^{2}y}{∂x∂p_{k}} )^{2} σ_{x}^{2}σ_{p_{k}}^{2} + ∑_{k} ∑_{l>k} (\frac{∂^{2}y}{∂p_{k}∂p_{l}} )^{2} σ_{p_{k}}^{2}σ_{p_{l}}^{2} +

+ \frac{∂y}{∂x} \frac{∂^{2}y}{∂x^{2}} E(ε_{x}^{3}) + ∑_{k} \frac{∂y}{∂p_{k}} \frac{∂^{2}y}{∂p_{k}^{2}} E(ε_{p_{k}}^{3}) .

The last relation can substantially be simplified using the approximate formulae E(ε_{x}^{4}) ≐ σ_{x}^{4}, E(ε_{p_{k}}^{4}) ≐ σ_{p_{k}}^{4} and E(ε_{x}^{3}) ≐ E(ε_{p_{k}}^{3}) ≐ 0, i.e.

D(ε) ≐ (\frac{∂y}{∂x})^{2} σ_{x}^{2} + σ_{y}^{2} + ∑_{k} (\frac{∂y}{∂p_{k}} )^{2} σ_{p_{k}}^{2} + ∑_{k} (\frac{∂^{2}y}{∂x∂p_{k}} )^{2} σ_{x}^{2}σ_{p_{k}}^{2} + ∑_{k} ∑_{l>k} (\frac{∂^{2}y}{∂p_{k}∂p_{l}} )^{2} σ_{p_{k}}^{2}σ_{p_{l}}^{2} ,

D(ε) ≐ V .

Unfortunately, both the relation for E(ε) and the relation for D(ε) are too complicated for most practical purposes anyway. The elements of the vector μ are

μ_{i} ≐ E(ε_{i}) ,

μ_{i} = \frac{1}{2} [(\frac{∂^{2}y}{∂x^{2}} )_{x=x_{i},y=y_{i}} σ_{x_{i}}^{2} + ∑_{k} (\frac{∂^{2}y}{∂p_{k}^{2}} )_{y=y_{i},p_{k}=p_{k}_{i}} σ_{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¶. If the quadratic instead of the linear Taylor approximation was applied, an estimate of the vector μ would have to be subtracted from the random vector y first.

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. The unknown linear and nonlinear parameters (variables of the vector of residuals f) are separated^{8}. The unknown nonlinear parameters are computed using the approximate Jacobi matrix^{9}

J_{2} = (1_{m} − AA^{+})J_{n−k}

(instead of the full-size Jacobi matrix J ≡ ( A J_{n−k}) of the function f with respect to all unknown parameters†), m is the number of points, n is the number of unknown parameters, k is the number of unknown linear parameters, A ∼ (m, k) and J_{n−k} ∼ (m, n − k) are the “linear” and “nonlinear” (left and right) submatrices of the matrix J ∼ (m, n)§.

When selected elements of the vector ϑ are the only unknown parameters, then a calculation is made of their unbiased estimates by the iterative 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.

The method used for optimization of the unknown parameters is based on an algorithm similar to the MDLS (\underline{m}odified \underline{d}amped \underline{l}east \underline{s}quares) method of Meyer and Roth^{10} with the partially restricted shift vector δ = −α(J^{T}J + λΩ)^{−1}J^{T}f (using the QR-decomposition of the Jacobi matrix‡), Ω is a diagonal matrix with the same diagonal elements as the matrix J^{T}J¶. Deliberate selection of a method of calculation of the unknown free concentrations of reactants of appropriate (i.e. not always the highest available) robustness considerably increases the efficiency of the method of optimization of the unknown parameters. Failure of the selected method (starting with the second iteration, the values computed in the previous iteration are taken for initial estimates of the unknown free concentrations) is evaluated as a step to a higher value on the hypersurface of sum of squares of residuals.

Moreover, several empirical tricks have been employed which significantly improve the “submethod” of optimization of the unknown stability constants. Probably the most effective of them is the change of all the negative unknown stability constants to positive ones (simply by their multiplication with −1) as soon as the λ modifier of step length and direction exceeds its upper limit. Not only that is why the OPIUM program rarely fails. Usually it finds the correct solution even from initial estimates of the unknown stability constants lying far beyond the limit of their physical sense.

¶ 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 jackknife and likelihood ratio analysis (interval estimate of mean, significance test). At each point, several 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 hypotheses of homogeneity, randomness, independence and normality of the standardized residuals.

† the efficiency of the separation of the unknown linear and nonlinear parameters compared with the solution of the full-size problem depends on the calibration function(s) used, initial estimates of the unknown nonlinear parameters and fitted data, but the computed results are always \underline{exactly the same}!

§ Strictly speaking, the matrix J_{n−k} is a function (of the unknown nonlinear parameters) composed of the “nonlinear” (right) submatrix of the matrix J (which is a function of both the unknown linear and the unknown nonlinear parameters) and the closed-form dependence of the unknown linear parameters on the unknown nonlinear parameters. Hence, optimum values of all the unknown linear parameters have to be computed (by the linear least squares method) in each iteration of the nonlinear least squares method (for actual values of the unknown nonlinear parameters and variance-covariance matrix V) before the matrix J_{n−k} is calculated and transformed to the matrix J_{2} by its left multiplication with the square matrix 1_{m} − AA^{+}.

‡ taking full advantage of the contingent block structure of its “linear” submatrix ¶ the λ modifier of step length and direction is initialized with the value 0.01

Four-level termination criterion has been adopted into the non-linear least squares method (only the 1.a level is applied in the case of iterative linear least squares method, “b” levels are used instead of “a” levels for testing of unknown free concentrations of reactants and stability constants computed in logarithmic scale):

1. a) |δ_{i}|/|p_{i}| ≤ ε_{1}, i = 1, …, n, b) |δ_{i}| ≤ − log(1 − ε_{1}), i = 1, …, n,

2. a) |g_{i}p_{i}| ≤ ε_{2}, i = 1, …, n, b) |g_{i}| ≤ ε_{2} ln 10, i = 1, …, n,

3. |cos α_{i}| = \frac{1}{2}|g_{i}|/\sqrt{j_{i}^{T}j_{i}f^{T}f} ≤ ε_{3}, i = 1, …, n,

4. the Jacobi matrix J has full rank, r(J) = n, n is the number of unknown parameters, ε_{1}, ε_{2}, ε_{3} are small real numbers, p is the vector of unknown parameters, δ is the vector of parameter shifts, f is the vector of residuals, j_{i} is the i-th column of the Jacobi matrix J of the function f (with respect to the unknown parameters), g = 2J^{T}f is the gradient of the function minimized ∥f∥^{2} = f^{T}f and α_{i} is the angle between j_{i} and f. This test ensures that the OPIUM program never terminates the nonlinear lest squares method computation at a false minimum or a saddle point.

Principal component analysis

The dependence of the measured quantity 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 .

¶ 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.

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 ,

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.

�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 defi-

nitely 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} .

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^{12}-Rissanen^{13}, Hannan-Quinn^{14} and Anděl-Perez-Negrao^{15} statistic) and stores it into a file§.

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

¶ To estimate the “effective” rank k of the matrix A (both the equations would hold exactly only if there were no measuring errors), sequence of residual standard deviations is calculated^{11} Then the 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 δ

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

1. Kvala M., Luke I., CHEMOMETRICS ’95, Pardubice, Czech Republic, 1995, Book of Abstracts p. 63 2. Havel J., Meloun M., Talanta 33 (1986) 525

3. Perrin D. D., Sayce I. G., Talanta 14 (1967) 833

4. Ginzburg G., Talanta 23 (1976) 149

5. Tripathi V. S., Talanta 33 (1986) 1015

6. De Robertis A., De Stefano C., Sammartano S., Anal. Chim. Acta 191 (1986) 1015

7. De Stefano C., Princi P., Rigano C., Sammartano S., Computers Chem. 13 (1989) 343

8. Golub G. H., Pereyra V., SIAM J. Numer. Anal. 10 (1973) 413

9. Maeder M., Zuberbühler A. D., Anal. Chem. 62 (1990) 2220

10. Meyer R. R., Roth D. M., J. Inst. Math. Appl. 9 (1972) 218

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

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

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

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

15. Andl J., Perez M. G., Negrao A. I., Kybernetika 17 (1981) 514