View source for Determination of equilibrium constants

[[Equilibrium constant]]s are determined in order to quantify [[chemical equilibria]]. When an equilibrium constant {{mvar|K}} is expressed as a concentration quotient,  
:<math>K=\frac{\mathrm{[S]} ^\sigma \mathrm{[T]}^\tau \cdots } {\mathrm{[A]}^\alpha \mathrm{[B]}^\beta \cdots }</math>
it is implied that the [[Activity (chemistry)|activity]] quotient is constant. For this assumption to be valid, equilibrium constants must be determined in a medium of relatively high [[ionic strength]]. Where this is not possible, consideration should be given to possible activity variation.
   
The equilibrium expression above is a function of the [[concentration]]s [A], [B] etc. of the chemical species in equilibrium. The equilibrium constant value can be determined if any one of these concentrations can be  measured. The general procedure is that the concentration in question is measured for a series of solutions with known [[analytical concentration]]s of the reactants. Typically, a [[titration]] is performed with one or more reactants in the titration vessel and one or more reactants in the burette. Knowing the analytical concentrations of reactants initially in the reaction vessel and in the burette, all analytical concentrations can be derived as a function of the volume (or mass) of titrant added.

The equilibrium constants may be derived by best-fitting of the experimental data with a chemical model of the equilibrium system.

== Experimental methods ==
There are four main experimental methods. For less commonly used methods, see Rossotti and Rossotti.<ref>{{cite book|first1=F.&nbsp;J.&nbsp;C.|last1=Rossotti|first2=H.|last2=Rossotti|title=The Determination of Stability Constants|publisher=McGraw-Hill|date=1961}}</ref> In all cases the range can be extended by using the [[#Competition method|competition method]]. An example of the application of this method can be found in [[palladium(II) cyanide]].

===  Potentiometric measurements ===
A free concentration [A] or activity {A} of a species A is measured by means of an [[ion selective electrode]] such as the [[glass electrode]]. If the electrode is calibrated using activity standards it is assumed that the [[Nernst equation]] applies in the form
:<math> E=E^0+\frac{RT}{nF}\ln\mathrm{\{A\}}</math>
where {{math|''E''<sup>0</sup>}} is the [[standard electrode potential]]. When [[Buffer solution|buffer]] solutions of known pH are used for calibration the [[pH meter|meter]] reading will be a pH.
:<math>\mathrm{pH}=\frac{nF}{RT}\left(E^0-E\right)</math>
At 298&nbsp;K, 1&nbsp;pH unit is approximately equal to 59&nbsp;mV.<ref>{{cite journal|url=http://www.iupac.org/publications/pac/1985/pdf/5703x0531.pdf|title=Definitions of pH scales, standard reference values, measurement of pH, and related terminology|journal=Pure Appl. Chem.|volume=57|pages=531–542|date=1985|doi=10.1351/pac198557030531|s2cid=14182410}}</ref>

When the electrode is calibrated with solutions of known concentration, by means of a strong acid–strong base titration, for example, a modified Nernst equation is assumed.
:<math>E=E^0 + s\log_{10}\mathrm{[A]}</math>
where {{mvar|s}} is an empirical slope factor. A solution of known hydrogen ion concentration may be prepared by standardization of a [[strong acid]] against [[borax]]. [[Azeotrope|Constant-boiling]] [[hydrochloric acid]] may also be used as a primary standard for hydrogen ion concentration.

==== Range and limitations ====
The most widely used electrode is the glass electrode, which is selective for the hydrogen ion. This is suitable for all [[Acid–base equilibrium|acid–base equilibria]]. {{math|log<sub>10</sub> ''β''}} values between about 2 and 11 can be measured directly by potentiometric [[titration]] using a [[glass electrode]]. This enormous range of stability constant values (ca. 100 to 10<SUP>11</sup>) is possible because of the logarithmic response of the electrode. The limitations arise because the Nernst equation breaks down at very low or very high pH.

When a glass electrode is used to obtain the measurements on which the calculated equilibrium constants depend, the precision of the calculated parameters is limited by secondary effects such as variation of [[liquid junction potential]]s in the electrode. In practice it is virtually impossible to obtain a precision for log β better than ±0.001.

=== Spectrophotometric measurements ===

==== Absorbance ====
It is assumed that the [[Beer–Lambert law]] applies.

:<math>A=l \sum {\varepsilon c}</math>

where {{mvar|l}} is the optical path length, {{mvar|ε}} is a molar absorbance at unit path length and {{mvar|c}} is a concentration. More than one of the species may contribute to the absorbance. In principle absorbance may be measured at one wavelength only, but in present-day practice it is common to record complete spectra.
===== Range and limitations =====
An upper limit on {{math|log<sub>10</sub> ''β''}} of 4 is usually quoted, corresponding to the precision of the measurements, but it also depends on how intense the effect is. Spectra of contributing species should be clearly distinct from each other

==== Fluorescence (luminescence) intensity ====
It is assumed that the scattered light intensity is a linear function of species’ concentrations.

:<math>I=\sum \varphi c </math>

where {{mvar|φ}} is a proportionality constant.
===== Range and limitations =====
 
The magnitude of the constant {{mvar|φ}} may be higher than the value of the molar extinction coefficient, ε, for a species. When this is so, the detection limit for that species will be lower. At high solute concentrations, fluorescence intensity becomes non-linear with respect to concentration due to self-absorption of the scattered radiation.

=== NMR chemical shift measurements ===
Chemical exchange is assumed to be rapid on the NMR time-scale. An individual chemical shift {{mvar|{{overbar|δ}}}} is the mole-fraction-weighted average of the shifts {{mvar|δ}} of nuclei in contributing species.
:<math>\bar {\delta} =\frac{\sum x_i \delta_i}{\sum x_i}</math>
Example: the p''K''<sub>a</sub> of the [[hydroxyl]] group in [[citric acid]] has been determined from [[Carbon-13 NMR|<sup>13</sup>C chemical shift]] data to be 14.4. Neither [[Ion-selective electrode|potentiometry]] nor [[ultraviolet–visible spectroscopy]] could be used for this determination.<ref>{{cite journal
| first1= Andre M. N.| last1= Silva
| first2= Xiaole| last2= Kong
| first3 =Robert C.| last3=Hider
| title= Determination of the p''K''<sub>a</sub> value of the hydroxyl group in the α-hydroxycarboxylates citrate, malate and lactate by <sup>13</sup>C NMR: implications for metal coordination in biological systems
| journal= Biometals| year = 2009| volume= 22| issue= 5| pages= 771–778| doi= 10.1007/s10534-009-9224-5| pmid=19288211
| s2cid= 11615864
}}</ref>

==== Range and limitations ====
Limited precision of chemical shift measurements also puts an upper limit of about 4 on {{math|log<sub>10</sub> ''β''}}. Limited to diamagnetic systems. [[Proton nuclear magnetic resonance|<sup>1</sup>H&nbsp;NMR]] cannot be used with solutions of compounds in <sup>1</sup>H<sub>2</sub>O.

=== Calorimetric measurements ===
Simultaneous measurement of {{mvar|K}} and {{math|[[Enthalpy|Δ''H'']]}} for 1:1 adducts is routinely carried out using [[isothermal titration calorimetry]]. Extension to more complex systems is limited by the availability of suitable software.

==== Range and limitations ====
Insufficient evidence is currently available.

=== The competition method ===
The competition method may be used when a stability constant value is too large to be determined by a direct method. It was first used by [[Gerold Schwarzenbach|Schwarzenbach]] in the determination of the stability constants of complexes of [[EDTA]] with metal ions.

For simplicity consider the determination of the stability constant <math>K_{AB}</math> of a binary complex, ''AB'', of a reagent ''A'' with another reagent ''B''.
:<math>K_{AB}=\frac{[AB]}{[A][B]}</math>    
where the [X] represents the concentration, at equilibrium, of a species X in a solution of given composition.

A ligand ''C'' is chosen which forms a weaker complex with ''A'' The stability constant, K<sub>AC</sub>, is small enough to be determined by a direct method. For example, in the case of EDTA complexes ''A'' is a  metal ion  and ''C'' may be a [[polyamine]] such as [[diethylenetriamine]].
:<math>K_{AC}=\frac{[AC]}{[A][C]}</math>    
The stability constant, ''K'' for the competition reaction 
:<math>AC + B \leftrightharpoons AB +C</math>
can be expressed as
:<math>K=\frac{[AB][C]}{[AC][B]}</math>
It follows that
: <math>K_{AB}=K \times K_{AC}</math>
where K is the stability constant for the competition reaction. Thus, the value of the stability constant <math>K_{AB}</math> may be derived from the experimentally determined values of ''K'' and <math>K_{AC}</math> .

== Computational methods ==
It is assumed that the collected experimental data comprise a set of data points. At each {{mvar|i}}th data point, the analytical concentrations of the reactants, {{math|''T''<sub>A</sub>(''i'')}}, {{math|''T''<sub>B</sub>(''i'')}} etc. are known along with a measured quantity, {{mvar|y<sub>i</sub>}}, that depends on one or more of these analytical concentrations. A general computational procedure has four main components:
# Definition of a chemical model of the equilibria
# Calculation of the concentrations of all the chemical species in each solution
# Refinement of the equilibrium constants
# Model selection
The value of the equilibrium constant for the formation of a 1:1 complex, such as a host-guest species, may be calculated with a dedicated spreadsheet application, Bindfit:<ref name=Hibbert>{{cite journal|first1=D.B.|last1=Hibbert|first2=P.|last2=Thordarson|journal=Chemical Communications|volume=52|issue=87|date=2017|doi=10.1039/c6cc03888c|title=The death of the Job plot, transparency, open science and online tools, uncertainty estimation methods and other developments in supramolecular chemistry data analysis|pages=12792–12805|pmid=27779264|doi-access=free}}</ref> In this case step 2 can be performed with a non-iterative procedure and the pre-programmed routine [[Solver]] can be used for step 3.

=== The chemical model ===
The chemical model consists of a set of chemical species present in solution, both the '''reactants''' added to the reaction mixture and the '''complex species''' formed from them. Denoting the reactants by A, B..., each ''complex species'' is specified by the [[stoichiometric coefficient]]s that relate the particular combination of ''reactants'' forming them.
:<chem>{\mathit p A} + \mathit q B \cdots <=> A_\mathit{p}B_\mathit{q} \cdots</chem>: <math chem>\beta_{pq\cdots}=\frac {[\ce{A}_p\ce{B}_q \cdots ]} {[\ce A]^p[\ce B]^q \cdots }</math>

When using general-purpose computer programs, it is usual to use [[Equilibrium constant#Cumulative and stepwise formation constants|cumulative]] [[Equilibrium constant#Association and dissociation constants|association constants]], as shown above. Electrical charges are not shown in general expressions such as this and are often omitted from specific expressions, for simplicity of notation. In fact, electrical charges have no bearing on the equilibrium processes other that there being a requirement for overall electrical neutrality in all systems.

With aqueous solutions the concentrations of proton (hydronium ion) and hydroxide ion are constrained by the self-dissociation of water.
:<chem>H2O <=> H+ + OH-</chem>: <math>K_\mathrm{W}^' =  \frac{[H^+][OH^-]}{[H_2O]} </math>

With dilute solutions the concentration of water is assumed constant, so the equilibrium expression is written in the form of the [[Self-ionization of water|ionic product]] of water.

:<math chem>K_\mathrm{W}=\ce{[H+]}[\ce{OH-}]\,</math>

When both H<sup>+</sup> and OH<sup>−</sup> must be considered as reactants, one of them is eliminated from the model by specifying that its concentration be derived from the concentration of the other. Usually the concentration of the hydroxide ion is given by

:<math chem>[\ce{OH-}]=\frac{K_\ce{W}}{[\ce{H+}]}\,</math>

In this case the equilibrium constant for the formation of hydroxide has the stoichiometric coefficients −1 in regard to the proton and zero for the other reactants. This has important implications for all protonation equilibria in aqueous solution and for [[Equilibrium constant#Hydrolysis constants|hydrolysis constants]] in particular.

It is quite usual to omit from the model those species whose concentrations are considered negligible. For example, it is usually assumed then there is no interaction between the reactants and/or complexes and the electrolyte used to maintain constant ionic strength or the buffer used to maintain constant pH. These assumptions may or may not be justified. Also, it is implicitly assumed that there are no other complex species present. When complexes are wrongly ignored a [[systematic error]] is introduced into the calculations.

Equilibrium constant values are usually estimated initially by reference to [[equilibrium constant#Data sources|data sources]].

=== Speciation calculations ===
A speciation calculation is one in which concentrations of all the species in an equilibrium system are calculated, knowing the analytical concentrations, T<sub>A</sub>, T<sub>B</sub> etc. of the reactants A, B etc. This means solving a set of nonlinear equations of mass-balance

:<math chem="">
\begin{align}
\ce{T_A} & = [\ce A]+\sum_{1,nk}p\beta_{pq \cdots}[\ce A]^p[\ce B]^q \cdots \\
\ce{T_B} & = [\ce B]+\sum_{1,nk}q\beta_{pq \cdots}[\ce A]^p[\ce B]^q \cdots \\
etc.
\end{align}
</math>

for the free concentrations [A], [B] etc. When the pH (or equivalent e.m.f., E).is measured, the free concentration of hydrogen ions, [H], is obtained from the measured value as  <blockquote><math>[\mathrm H]=10^{-\mathrm{pH}}</math> or <math>[\mathrm H]=e^\mathrm{{ -\frac{nF}{RT}}(E-E^0)  }</math></blockquote>and only the free concentrations of the other reactants are calculated. The concentrations of the complexes are derived from the free concentrations via the chemical model.

Some authors<ref>{{cite journal|first1=R.&nbsp;J.|last1=Motekaitis|first2=A.&nbsp;E.|last2=Martell|journal=Can. J. Chem.|date=1982|volume=60|pages=2403–2409|doi=10.1139/v82-347|title=BEST — A new program for rigorous calculation of equilibrium parameters of complex multicomponent systems|issue=19}}</ref><ref name=pgp1990a>{{cite journal|first1=P.&nbsp;G.|last1=Potvin|journal=Can. J. Chem.|date=1990|volume=68|pages=2198–2207|doi=10.1139/v90-337|title=Modelling complex solution equilibria. I. Fast, worry-free least-squares refinement of equilibrium constants|issue=12|doi-access=free}}</ref> include the free reactant terms in the sums by declaring ''identity'' (unit) {{mvar|β}} constants for which the stoichiometric coefficients are 1 for the reactant concerned and zero for all other reactants. For example, with 2 reagents, the mass-balance equations assume the simpler form.

:<math chem>
\begin{align}
T_\ce{A} & = \sum_{0,nk}p\beta_{pq}[\ce A]^p[\ce B]^q \\[4pt]
T_\ce{B} & = \sum_{0,nk}q\beta_{pq}[\ce A]^p[\ce B]^q \\
\end{align}
</math>
:<math chem="">\beta_{10}= \beta_{01} = 1</math>
In this manner, all chemical species, ''including the free reactants'', are treated in the same way, having been ''formed'' from the combination of reactants that is specified by the stoichiometric coefficients.

In a titration system the analytical concentrations of the reactants at each titration point are obtained from the initial conditions, the burette concentrations and volumes. The analytical (total) concentration of a reactant R at the {{mvar|i}}th titration point is given by
:<math chem>T_\ce{R}=\frac{\ce{R}_0+v_i\ce{[R]}}{v_0+v_i}</math>

where R<sub>0</sub> is the initial '''amount''' of R in the titration vessel, {{math|''v''<sub>0</sub>}} is the initial volume, [R] is the '''concentration''' of R in the burette and {{mvar|v<sub>i</sub>}} is the volume added. The burette concentration of a reactant not present in the burette is taken to be zero.

In general, solving these nonlinear equations presents a formidable challenge because of the huge range over which the free concentrations may vary. At the beginning, values for the free concentrations must be estimated. Then, these values are refined, usually by means of [[Newton's method|Newton–Raphson]] iterations. The logarithms of the free concentrations may be refined rather than the free concentrations themselves. Refinement of the logarithms of the free concentrations has the added advantage of automatically imposing a non-negativity constraint on the free concentrations. Once the free reactant concentrations have been calculated, the concentrations of the complexes are derived from them and the equilibrium constants.

Note that the free reactant concentrations can be regarded as implicit parameters in the equilibrium constant refinement process. In that context the values of the free concentrations are constrained by forcing the conditions of mass-balance to apply at all stages of the process.

=== Equilibrium constant refinement ===
The objective of the refinement process is to find equilibrium constant values that give the best fit to the experimental data. This is usually achieved by minimising an [[objective function]], {{mvar|U}}, by the method of non-linear [[least-squares]]. First the residuals are defined as

:<math chem>r_i=y_i^\text{obs}-y_i^\text{calc}</math>

Then the most general objective function is given by

:<math chem>U=\sum_i\sum_j r_i W_{ij} r_j\,</math>

The matrix of weights, {{math|'''W'''}}, should be, ideally, the inverse of the [[Covariance matrix|variance-covariance matrix]] of the observations. It is rare for this to be known. However, when it is, the [[expectation value]] of U is one, which means that the data are fitted ''within experimental error''. Most often only the diagonal elements are known, in which case the objective function simplifies to

:<math chem>U=\sum_i W_{ii}r_i^2</math>

with {{math|''W<sub>ij</sub>'' {{=}} 0}} when {{math|''j'' ≠ ''i''}}. Unit weights, {{math|''W<sub>ii</sub>'' {{=}} 1}}, are often used but, in that case, the expectation value of {{math|U}} is the [[root mean square]] of the experimental errors.

The minimization may be performed using the [[Gauss–Newton method]]. Firstly the objective function is linearised by approximating it as a first-order [[Taylor series]] expansion about an initial parameter set, {{math|'''p'''}}.
:<math chem>U=U^0+\sum_i \frac{\partial U}{\partial p_i}\delta p_i</math>
The increments {{math|δ''p<sub>i</sub>''}} are added to the corresponding initial parameters such that {{mvar|U}} is less than {{math|''U''<sup>0</sup>}}. At the minimum the derivatives {{math|{{sfrac|∂''U''|∂''p<sub>i</sub>''}}}}, which are simply related to the elements of the [[Jacobian matrix and determinant|Jacobian]] matrix, {{math|'''J'''}}

:<math>J_{jk}=\frac{\partial y_j^\mathrm{calc}}{\partial p_k}</math>

where {{mvar|p<sub>k</sub>}} is the {{mvar|k}}th parameter of the refinement, are equal to zero. One or more equilibrium constants may be parameters of the refinement. However, the measured quantities (see above) represented by {{math|'''y'''}} are not expressed in terms of the equilibrium constants, but in terms of the species concentrations, which are implicit functions of these parameters. Therefore, the Jacobian elements must be obtained using [[implicit differentiation]].

The parameter increments {{math|δ'''p'''}} are calculated by solving the [[normal equations]], derived from the conditions that {{math|{{sfrac|∂'''U'''|∂'''p'''}} {{=}} 0}} at the minimum.
:<math chem>{ \left(J^\mathrm{T} W J\right) \delta p=J^\mathrm{T} W r }</math>

The increments {{math|δ'''p'''}} are added iteratively to the parameters
:<math chem>\mathbf{p}^{n+1}=\mathbf{p}^n +\delta \mathbf{p}</math>
where {{mvar|n}} is an iteration number. The species concentrations and {{math|''y''<sup>calc</sup>}} values are recalculated at every data point. The iterations are continued until no significant reduction in {{mvar|U}} is achieved, that is, until a convergence criterion is satisfied. If, however, the updated parameters do not result in a decrease of the objective function, that is, if divergence occurs, the increment calculation must be modified. The simplest modification is to use a fraction, {{mvar|f}}, of calculated increment, so-called shift-cutting.
:<math chem>\mathbf{p}^{n+1}=\mathbf{p}^n +f \mathbf{\delta p}</math>
In this case, the direction of the shift vector, {{math|δ'''p'''}}, is unchanged. With the more powerful [[Levenberg–Marquardt algorithm]], on the other hand, the shift vector is rotated towards the direction of [[steepest descent]], by modifying the normal equations, 
:<math chem>\mathbf{ \left(J^\mathrm{T} W J +\lambda I\right)\delta p=J^\mathrm{T} W r }</math>
where {{mvar|λ}} is the Marquardt parameter and {{math|'''I'''}} is an identity matrix. Other methods of handling divergence have been proposed.<ref name="pgp1990a" />

A particular issue arises with NMR and spectrophotometric data. For the latter, the observed quantity is absorbance, {{mvar|A}}, and the Beer–Lambert law can be written as
:<math chem>A_\lambda=l\sum(\varepsilon_{pq..})_\lambda c_{pq..}</math>

It can be seen that, assuming that the concentrations, c, are known, that absorbance, {{mvar|A}}, at a given wavelength, <math chem>\lambda</math>, and path length <math chem>l</math>, is a linear function of the molar absorptivities, {{mvar|ε}}. With 1&nbsp;cm path-length, in matrix notation

:<math chem>\mathbf{A}=\boldsymbol{\varepsilon} \mathbf{C} \, </math>

There are two approaches to the calculation of the unknown molar absorptivities
:(1) The {{math|'''ε'''}} values are considered parameters of the minimization and the Jacobian is constructed on that basis. However, the {{math|'''ε'''}} values themselves are calculated at each step of the refinement by linear least-squares:
::<math chem>\boldsymbol{\varepsilon} = \mathbf{\left(C^\mathrm{T}C\right)^{-1}C^\mathrm{T}A }</math>
:using the refined values of the equilibrium constants to obtain the speciation. The matrix
::<math chem>\mathbf{\left(C^TC\right)^{-1}C^T}</math>
:is an example of a [[pseudo-inverse]].
:Golub and Pereyra<ref name="Golub">{{cite journal|first1=G.&nbsp;H.|last1=Golub|first2=V.|last2=Pereyra|journal=SIAM J. Numer. Anal.|volume=10|pages=413–432|date=1973|doi=10.1137/0710036|title=The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems Whose Variables Separate|issue=2|bibcode=1973SJNA...10..413G}}</ref> showed how the pseudo-inverse can be differentiated so that parameter increments for both molar absorptivities and equilibrium constants can be calculated by solving the normal equations.
:(2) The Beer–Lambert law is written as
::<math chem>\mathbf{\boldsymbol\varepsilon}_\lambda= \mathbf{A}^{-1}_\lambda \mathbf{C} \, </math>
:The unknown molar absorbances of all "coloured" species  are found by using the non-iterative method of linear least-squares, one wavelength at a time. The calculations are performed once every refinement cycle, using the stability constant values obtaining at that refinement cycle to calculate species' concentration values in the matrix <math chem> \mathbf{C} </math>.

==== Parameter errors and correlation ====
In the region close to the minimum of the objective function, {{mvar|U}}, the system approximates to a linear least-squares system, for which
:<math chem>\mathbf{p=\left(J^\mathrm{T}WJ\right)^{-1}J^\mathrm{T}Wy^\mathrm{obs}}</math>
Therefore, the parameter values are (approximately) linear combinations of the observed data values and the errors on the parameters, {{math|'''p'''}}, can be obtained by [[error propagation]] from the observations, {{math|'''y'''<sup>obs</sup>}}, using the linear formula. Let the [[variance-covariance matrix]] for the observations be denoted by {{math|'''Σ<sup>y</sup>'''}} and that of the parameters by {{math|'''Σ<sup>p</sup>'''}}. Then,
:<math chem>\mathbf{\Sigma^p=\left(J^\mathrm{T}WJ\right)^{-1}J^\mathrm{T}W \Sigma^y W^\mathrm{T}J(J^\mathrm{T}WJ)^{-1}}</math>
When {{math|'''W''' {{=}} ('''Σ<sup>y</sup>''')<sup>−1</sup>}}, this simplifies to
:<math chem>\mathbf{\Sigma^p=\left(J^\mathrm{T}WJ\right)^{-1}}</math>
In most cases the errors on the observations are un-correlated, so that {{math|'''Σ<sup>y</sup>'''}} is [[diagonal matrix|diagonal]].
If so, each weight should be the [[Multiplicative inverse|reciprocal]] of the variance of the corresponding observation. For example, in a [[potentiometric titration]], the weight at a titration point, {{mvar|k}}, can be given by
:<math chem>W_k= \frac{1}{\sigma^2_E+\left( \frac{\partial E}{\partial v} \right)^2_k\sigma^2_v} </math>
where {{mvar|σ<sub>E</sub>}} is the error in electrode potential or pH, {{math|<big><big>(</big></big>{{sfrac|∂''E''|∂''v''}}<big><big>)</big></big>{{su|b=''k''}}}} is the slope of the [[titration curve]] and {{mvar|σ<sub>v</sub>}} is the error on added volume.

When unit weights are used ({{math|'''W''' {{=}} '''I'''}}, {{math|'''p''' {{=}} ('''J'''<sup>T</sup>'''J''')<sup>−1</sup>'''J'''<sup>T</sup>'''y'''}}) it is implied that the experimental errors are uncorrelated and all equal: {{math|'''Σ<sup>y</sup>''' {{=}} ''σ''<sup>2</sup>'''I'''}}, where {{math|''σ''<sup>2</sup>}} is known as the variance of an observation of unit weight, and {{math|'''I'''}} is an [[identity matrix]]. In this case {{math|''σ''<sup>2</sup>}} is approximated by
:<math chem>\sigma^2 = \frac{U}{n_\mathrm{d}-n_\mathrm{p}}</math>
where {{mvar|U}} is the minimum value of the objective function and {{math|''n''<sub>d</sub>}} and {{math|''n''<sub>p</sub>}} are the number of data and parameters, respectively.
:<math chem>\mathbf{\Sigma^p}=\frac{U}{n_\mathrm{d}-n_\mathrm{p}}\left(\mathbf{J}^\mathrm{T}\mathbf{J}\right)^{-1}</math>
In all cases, the [[variance]] of the parameter {{mvar|p<sub>i</sub>}} is given by {{mvar|Σ{{su|p=p|b=ii}}}} and the covariance between parameters {{mvar|p<sub>i</sub>}} and {{mvar|p<sub>j</sub>}} is given by {{mvar|Σ{{su|p=p|b=ij}}}}. [[Standard deviation]] is the square root of variance. These error estimates reflect only [[random errors]] in the measurements. The true uncertainty in the parameters is larger due to the presence of [[systematic errors]]—which, by definition, cannot be quantified.

Note that even though the observations may be uncorrelated, the parameters are always [[Pearson product-moment correlation coefficient|correlated]].

==== Derived constants====
When cumulative constants have been refined it is often useful to derive stepwise constants from them. The general procedure is to write down the defining expressions for all the constants involved and then to equate concentrations. For example, suppose that one wishes to derive the pKa for removing one proton from a tribasic acid, LH<sub>3</sub>, such as [[citric acid]].
:<math chem>\begin{align}
\ce{L^3-}+ \ce{ H+ <=> }\ \ce{LH^2-} &:\ [\ce{LH^2-}] =\beta_{11} [\ce{L^3-}] [\ce{H+}]\\
\ce{L^3-}+ \ce{2H+ <=> }\ \ce{LH2^-} &:\ [\ce{LH2^-}] =\beta_{12} [\ce{L^3-}] [\ce{H+}]^2\\
\ce{L^3-}+ \ce{3H+ <=> }\ \ce{LH3}   &:\ [\ce{LH3}]   =\beta_{13} [\ce{L^3-}] [\ce{H+}]^3
\end{align}</math>
The stepwise ''association'' constant for formation of LH<sub>3</sub> is given by
:<math chem>\ce{{LH2^-} + H+ <=> LH3\ ; \quad\ [LH3]}=K[\ce{LH2^-}][\ce{H+}]</math>
Substitute the expressions for the concentrations of LH<sub>3</sub> and {{chem|LH|2|−}} into this equation
:<math chem>\beta_{13}[\ce{L^3-}][\ce{H+}]^3=K\beta_{12}[\ce{L^3-}][\ce{H+}]^2[\ce{H+}]</math>
whence
:<math>\beta_{13}=K\beta_{12}; K=\frac{\beta_{13}}{\beta_{12}} \,</math>
and since {{math|p''K''<sub>a</sub> {{=}} −log<sub>10</sub> {{sfrac|1|''K''}}}} its value is given by
:<math chem>\ce{p}K_\ce{a1} = \log_{10} \beta_{13}-\log_{10} \beta_{12}\, </math>
:<math chem>\ce{p}K_\ce{a2} = \log_{10} \beta_{12}-\log_{10} \beta_{11}\, </math>
:<math chem>\ce{p}K_\ce{a3} = \log_{10} \beta_{11}\, </math>
Note the reverse numbering for pK and log β. When calculating the error on the stepwise constant, the fact that the cumulative constants are [[correlated]] must accounted for. By [[error propagation]]
:<math>\sigma^2_K=\sigma^2_{\beta_{12}}+\sigma^2_{\beta_{13}}-2 \sigma_{\beta_{12}} \sigma_{\beta_{13}}\rho_{12,13}\,</math>
and
:<math>\sigma_{\log_{10} K}=\frac{\sigma_K}{K}</math>

=== Model selection ===
Once a refinement has been completed the results should be checked to verify that the chosen model is acceptable. generally speaking, a model is acceptable when the data are fitted within experimental error, but there is no single criterion to use to make the judgement. The following should be considered.

==== The objective function ====
When the weights have been correctly derived from estimates of experimental error, the [[expectation value]] of {{math|{{sfrac|''U''|''n''<sub>d</sub> − ''n''<sub>p</sub>}}}} is 1.<ref name=Hamilton>{{cite book|first=W.&nbsp;C.|last=Hamilton|title=Statistics in Physical Science|url=https://archive.org/details/statisticsinphys0000hami|url-access=registration|publisher=Ronald Press|location=New York, NY|date=1964}}</ref> It is therefore very useful to estimate experimental errors and derive some reasonable weights from them as this is an absolute indicator of the goodness of fit.

When unit weights are used, it is implied that all observations have the same variance. {{math|{{sfrac|''U''|''n''<sub>d</sub> − ''n''<sub>p</sub>}}}} is expected to be equal to that variance.

==== Parameter errors ====
One would want the errors on the stability constants to be roughly commensurate with experimental error. For example, with pH titration data, if pH is measured to 2 decimal places, the errors of {{math|log<sub>10</sub> ''β''}} should not be much larger than 0.01. In exploratory work where the nature of the species present is not known in advance, several different chemical models may be tested and compared. There will be models where the uncertainties in the best estimate of an equilibrium constant may be somewhat or even significantly larger than {{math|''σ''<sub>pH</sub>}}, especially with those constants governing the formation of comparatively minor species, but the decision as to how large is acceptable remains subjective. The decision process as to whether or not to include comparatively uncertain equilibria in a model, and for the comparison of competing models in general, can be made objective and has been outlined by Hamilton.<ref name=Hamilton />

==== Distribution of residuals ====
[[Image:NiGly titration.jpg|right]]
At the minimum in {{mvar|U}} the system can be approximated to a linear one, the residuals in the case of unit weights are related to the observations by
:<math chem>\mathbf{r=y^\mathrm{obs}-J \left(J^\mathrm{T}T \right)^{-1}J^\mathrm{T} y^\mathrm{obs}}</math>
The [[symmetric matrix|symmetric]], [[idempotent]] matrix {{math|'''J'''('''J'''<sup>T</sup>'''T''')<sup>−1</sup>'''J'''}} is known in the statistics literature as the [[hat matrix]], {{math|'''H'''}}. Thus,
:<math chem>\mathbf{r=\left(I-H \right) y^\mathrm{obs}}</math>
and
:<math chem>\mathbf{M^r=\left(I-H \right) M^y \left(I-H \right)}</math>
where {{math|'''I'''}} is an [[identity matrix]] and {{math|'''M<sup>r</sup>'''}} and {{math|'''M<sup>y</sup>'''}} are the variance-covariance matrices of the residuals and observations, respectively. This shows that even though the observations may be uncorrelated, the residuals are always correlated.

The diagram at the right shows the result of a refinement of the stability constants of [[nickel|Ni]](Gly)<sup>+</sup>, Ni(Gly)<sub>2</sub> and {{chem|Ni(Gly)|3|-}} (where GlyH = [[glycine]]). The observed values are shown a blue diamonds and the species concentrations, as a percentage of the total nickel, are superimposed. The residuals are shown in the lower box. The residuals are not distributed as randomly as would be expected. This is due to the variation of [[liquid junction potential]]s and other effects at the glass/liquid interfaces. Those effects are very slow compared to the rate at which equilibrium is established.

==== Physical constraints ====
Some physical constraints are usually incorporated in the calculations. For example, all the concentrations of free reactants and species must have positive values and association constants must have positive values.

With spectrophotometric data the calculated molar absorptivity (or emissivity) values should all be positive. Most computer programs do not impose this constraint on the calculations.

==== Chemical constraints ====
When determining the stability constants of metal-ligand complexes, it is common practice to fix ligand protonation constants at values that have been determined using data obtained from metal-free solutions. Hydrolysis constants of metal ions are usually fixed at values which were obtained using ligand-free solutions. When determining the stability constants for ternary complexes, M<sub>p</sub>A<sub>q</sub>B<sub>r</sub> it is common practice the fix the values for the corresponding binary complexes M<sub>p′</sub>A<sub>q′</sub> and M<sub>p′′</sub>B<sub>q′′</sub>, at values which have been determined in separate experiments. Use of such constraints reduces the number of parameters to be determined, but may result in the calculated errors on refined stability constant values being under-estimated.

==== Other models ====

If the model is not acceptable, a variety of other models should be examined to find one that best fits the experimental data, within experimental error. The main difficulty is with the so-called minor species. These are species whose concentration is so low that the effect on the measured quantity is at or below the level of error in the experimental measurement. The constant for a minor species may prove impossible to determine if there is no means to increase the concentration of the species.
.

==Thermodynamic principles of host–guest interactions==
The thermodynamics of the host- guest interaction can be assessed by NMR spectroscopy, UV/visible spectroscopy, and isothermal titration calorimetry.<ref>{{cite journal | last1 = Piñeiro | first1 = Á. | last2 = Banquy | first2 = X. | last3 = Pérez-Casas | first3 = S. | last4 = Tovar | first4 = É. | last5 = García | first5 = A. | last6 = Villa | first6 = A. | last7 = Amigo | first7 = A. | last8 = Mark | first8 = A. E. | last9 = Costas | first9 = M. | title = On the Characterization of Host–Guest Complexes: Surface Tension, Calorimetry, and Molecular Dynamics of Cyclodextrins with a Non-ionic Surfactant | journal = Journal of Physical Chemistry B | volume = 111 | year = 2007 | pages = 4383–92 | issue = 17| pmid = 17428087| doi = 10.1021/jp0688815}}</ref> Quantitative analysis of binding constant values provides useful thermodynamic information.{{Citation needed|date=January 2025}}<!-- undefined <ref name = textbook/> -->

An [[association constant]], <math>K^\ominus_a</math> can be defined by the expression
:<math>K^\ominus_a = \frac{\{HG\}}{\{H\}\{G\}} = \frac{[HG]}{[H][G]} \times \Gamma </math>
where {HG} is the [[thermodynamic activity]] of the complex at equilibrium. {H} represents the activity of the host and {G} the activity of the guest. The quantities <math>[HG]</math>, <math>[H]</math> and <math>[G]</math> are the corresponding concentrations and <math>\Gamma</math> is a quotient of [[activity coefficient]]s.

In practice the equilibrium constant is usually defined in terms of concentrations. 
:<math>K_a =\frac{[HG]}{[H][G]}</math>
When this definition is used, it is implied that the quotient of [[activity coefficient]]s has a numerical value of one. It then appears that the equilibrium constant, <math>K_A</math> has the [[Dimensional analysis|dimension]] 1/concentration, but that cannot be true since the standard Gibbs free energy change, <math>\Delta G^\ominus</math> is proportional to the logarithm of <math>K_A</math>.
:<math>\Delta G^\ominus = -RT \ln{K_A^\ominus} </math>
This apparent paradox is resolved when the dimension of <math>\Gamma</math> is ''defined'' to be the reciprocal of the dimension of the quotient of concentrations. The implication is that <math>\Gamma</math> is regarded as having a constant value under all relevant experimental conditions. Nevertheless it is common practice to attach a dimension, such as millimole per litre or micromole per litre, to a value of ''K'' that has been determined experimentally.

A Large <math>K_a</math> value indicates that host and guest molecules interact strongly to form the host–guest complex.
<!-- [[File:HGfreeenergydiagram2.jpg|thumb|250px|right|Host–guest complex formation. X axis = reaction coordinate. Y axis = free energy. ]] 
This equilibrium constant has the [[Dimensional analysis|dimension]] 1/concentration and so its logarithm does not exist. To overcome this problem each concentration term should be divided by a standard concentration. When the value of each standard concentration is 1, the value of <math>K_a</math> is unchanged. Put another way, when the equilibrium constant is calculated as a quotient of concentrations it is implicitly assumed that the [[activity coefficient]] of each chemical species has a value of 1, which is effectively a definition of the thermodynamic [[standard state]]s.-->

===Determination of binding constant values and kinetic constant===
{{main| Determination of equilibrium constants}}

====Simple host–guest complexation====
When the host and guest molecules combine to form a single complex, the equilibrium is represented as
:<math>H + G \leftrightharpoons HG</math>
and the equilibrium constant, K, is defined as
:<math>K=\frac{[HG]}{[H] [G]}</math>
where [X] denotes the concentration of a chemical species X (all activity coefficients are assumed to have a numerical values of 1).
The mass-balance equations, at any data point, 
:<math>T_H = [H] + K [H] [G]</math> 
:<math>T_G = [G] + K [H] [G]</math> 
where <math>T_G</math> and <math>T_H</math> represent the total concentrations, of host and guest, can be reduced to a single quadratic equation in, say, [G] and so can be solved analytically for any given value of K. The concentrations [H] and [HG] can then derived. 
:<math>[H] = T_H - T_G + [G]</math>
:<math>[HG] = K [H] [G]</math>
The next step in the calculation is to calculate the value, <math>X^{calc}_i</math>, of a quantity corresponding to the quantity observed <math>X^{obs}_i</math>. Then, a sum of squares, U, over all data points, np, can be defined as 
:<math>U=\sum_{i=1,np} (X^{obs}_i -X^{calc}_i)^2</math>
and this can be minimized with respect to the stability constant value, K, and a parameter such the chemical shift of the species HG (nmr data) or its molar absorbency (uv/vis data). The minimization can be performed in a spreadsheet application such as [[Microsoft Excel|EXCEL]] by using the in-built SOLVER utility.

This procedure is applicable to 1:1 adducts.

====General complexation reaction====
For each equilibrium involving a host, H, and a guest G
:<math>pH + qG \leftrightharpoons H_pG_q</math>
the equilibrium constant, <math>\beta _{pq}</math>, is defined as 
:<math>\beta _{pq}=\frac{[H_p G_q]}{[H]^p [G]^q}</math> 
The values of the free concentrations, <math>[H]</math> and <math>[G]</math> are obtained by solving the equations of mass balance with known or estimated values for the stability constants.
:<math>T_H = [H] + \sum p\beta_{pq}[H]^p [G]^q</math> 
:<math>T_G = [G] + \sum q\beta_{pq}[H]^p [G]^q</math> 
Then, the concentrations of each complex species may also be calculated as <math>[H_p G_q] =\beta_{pq}[H]^p [G]^q </math>. The relationship between a species' concentration and the measured quantity is specific for the measurement technique, as indicated in each section above. Using this relationship, the set of parameters, the stability constant values and values of properties such as molar absorptivity or specified chemical shifts, may be refined by a non-linear [[least-squares]] refinement process. For a more detailed exposition of the theory see [[Determination of equilibrium constants]]. 
Some dedicated computer programs are listed at [[Determination of equilibrium constants#Implementations|Implementations]].

===Cooperativity===
In cooperativity, the initial ligand binding affects the host's affinity for subsequent ligands. In positive cooperativity, the first binding event enhances the affinity of the host for another ligand. Examples of positive and negative cooperativity are [[hemoglobin]] and aspartate receptor, respectively.<ref>{{cite journal | title = The structural basis of negative cooperativity: receptors and enzymes| journal = [[Current Opinion in Structural Biology]]| year = 1996| volume = 6 | pages = 757–761| doi = 10.1016/S0959-440X(96)80004-2 | author1 = Koshland, D | pmid=8994875 | issue = 6}}</ref>
[[File:Cooperativity3.png|thumb|300px|right|General Host–Guest Binding. <span style="color:black">(1.)</span> Guest A binding <span style="color:black">(2.)</span> Guest B binding. <span style="color:black">(3.)</span> Positive Cooperativity Guest A–B binding. <span style="color:black">(4.)</span> Negative Cooperativity Guest A–B binding]]
The thermodynamic properties of cooperativity have been studied in order to define mathematical parameters that distinguish positive or negative cooperativity. The traditional Gibbs free energy equation states: <math>\Delta G = \Delta H - T \Delta S \ </math>. However, to quantify cooperativity in a host–guest system, the binding energy needs to be considered. The schematic on the right shows the binding of A, binding of B, positive cooperative binding of A–B, and lastly, negative cooperative binding of A–B. Therefore, an alternate form of the Gibbs free energy equation would be
:<math>\Delta G_S^\circ = \Delta G_A^\circ + \Delta G_B^\circ - \Delta G_{AB}^\circ </math>
:<math>\Delta H_S^\circ = \Delta H_A^\circ + \Delta H_B^\circ - \Delta H_{AB}^\circ </math>
:<math>\ T\Delta G_S^\circ = T\Delta H_A^\circ + T\Delta H_B^\circ - T\Delta S_{AB}^\circ </math>
where:
:<math>\Delta G_A^\circ </math> = free energy of binding A
:<math>\Delta G_B^\circ </math> = free energy of binding B
:<math>\Delta G_S^\circ </math> = free energy of binding for A and B tethered
:<math>\Delta G_{AB}^\circ </math> = sum of the free energies of binding

It is considered that if <math>\Delta G_S^\circ </math> more than the sum of <math>\Delta G_A^\circ </math> and <math>\Delta G_B^\circ </math>, it is positively cooperative. If <math>\Delta G_S^\circ </math> is less, then it is negatively cooperative.<ref>{{cite journal | doi = 10.1073/pnas.78.7.4046 | title = On the attribution and additivity of binding energies | journal = [[Proceedings of the National Academy of Sciences, USA]] | year = 1981 | volume = 78 | issue = 7 | pages = 4046–4050 | pmid = 16593049 | author1 = Jencks, W. P. | pmc = 319722 |bibcode = 1981PNAS...78.4046J| doi-access = free }}</ref>
Host–guest chemistry is not limited to receptor-lingand interactions. It is also demonstrated in ion-pairing systems. Such interactions are studied in an aqueous media utilizing synthetic organometallic hosts and organic guest molecules. For example, a poly-cationic receptor containing copper (the host) is coordinated with molecules such as tetracarboxylates, tricarballate, aspartate, and acetate (the guests). This study illustrates that [[entropy]] rather than [[enthalpy]] determines the binding energy of the system leading to negative cooperativity. The large change in entropy originates from the displacement of solvent molecules surrounding the ligand and the receptor. When multiple acetates bind to the receptor, it releases more water molecules to the environment than a tetracarboxylate. This led to a decrease in free energy implying that the system is cooperating negatively.<ref>{{cite journal | title = Studies into the Thermodynamic Origin of Negative Cooperativity in Ion-Pairing Molecular Recognition | journal = [[Journal of the American Chemical Society]] | year = 2003 | volume = 125 | pages = 10963–10970 | doi = 10.1021/ja030265o | author1 = Dobrzanska, L | author2 = Lloyd, G | author3 = Esterhuysen, C | author4 = Barbour, L | issue = 36| pmid = 12952478 }}</ref> In a similar study, utilizing guanidinium and Cu(II) and polycarboxylate guests, it is demonstrated that positive cooperatively is largely determined by enthalpy.<ref>{{cite journal | title = A cationic host displaying positive cooperativity in water | journal = [[Proceedings of the National Academy of Sciences, USA]] | year = 2007 | volume = 104 | pages = 6538–6543 | doi = 10.1073/pnas.0609144104 | author1 = Hughes, A. | author2 = Anslyn, E | pmid=17420472 | issue = 16 | pmc = 1871821|bibcode = 2007PNAS..104.6538H| doi-access = free }}</ref> In addition to thermodynamic studies, host–guest chemistry also has biological applications.

== Implementations ==
Some simple systems are amenable to spreadsheet calculations.<ref name=Hibbert/><ref>{{cite book |last=Billo |first=E. Joseph |title=Excel for Chemists: A Comprehensive Guide |edition=3rd |year=2011 |publisher=Wiley-VCH |isbn= 978-0-470-38123-6}}</ref>

A large number of general-purpose computer programs for equilibrium constant calculation have been published. See <ref name="HQ">{{cite journal|first1=P.|last1=Gans|first2=A.|last2=Sabatini|first3=A.|last3=Vacca|title=Investigation of equilibria in solution. Determination of equilibrium constants with the HYPERQUAD suite of programs|journal=Talanta|volume=43|pages=1739–1753|date=1996|issue=10|doi=10.1016/0039-9140(96)01958-3|pmid=18966661}}</ref> for a bibliography. The most frequently used programs are:
* Potentiometric data: [http://www.hyperquad.co.uk/hq2000.htm Hyperquad], BEST<ref>{{cite book|first1=A.&nbsp;E.|last1=Martell|first2=R.&nbsp;J.|last2=Motekaitis|title=The Determination and Use of Stability Constants|publisher=Wiley-VCH|date=1992|isbn=0471188174}}</ref> PSEQUAD,<ref name=Leggett>{{cite book|editor-first=D.&nbsp;J.|editor-last=Leggett|title=Computational Methods for the Determination of Formation Constants|publisher=Plenum Press|date=1985|isbn=978-0-306-41957-7}}</ref> [http://jplusconsulting.com/products/reactlab-ph/ ReactLab pH PRO]
* Spectrophotometric data:[http://www.hyperquad.co.uk/HypSpec.htm HypSpec], SQUAD,<ref name="Leggett"/> Specfit,<ref>{{cite journal|first1=H.|last1=Gampp|first2=M.|last2=Maeder|first3=C.&nbsp;J.|last3=Mayer|first4=A.|last4=Zuberbühler|journal=Talanta|volume=32|issue=95|date=1985|doi=10.1016/0039-9140(85)80035-7|title=Calculation of equilibrium constants from multiwavelength spectroscopic data—IMathematical considerations|pages=95–101|pmid=18963802}}</ref> [http://jplusconsulting.com/products/reactlab-equilibria/ ReactLab EQUILIBRIA] 
* NMR data [http://www.hyperquad.co.uk/hypnmr.htm HypNMR], [http://www.nuigalway.ie/chem/Mike/wineqnmr.htm EQNMR ] {{Webarchive|url=https://web.archive.org/web/20190714173527/http://www.nuigalway.ie/chem/Mike/wineqnmr.htm |date=2019-07-14 }}
*Calorimetric data [http://www.hyperquad.co.uk/hypdeltah.htm HypΔH]. [https://www.affinimeter.com/site/itc/ Affinimeter] <br>Commercial [[Isothermal titration calorimeter]]s are usually supplied with software with which an equilibrium constant and [[enthalpy|standard formation enthalpy]] for the formation of a 1:1 adduct can be obtained. Some software for handling more complex equilibria may also be supplied.

== References ==
{{reflist|30em}}{{Chemical equilibria}}
{{DEFAULTSORT:Determination Of Equilibrium Constants}}
[[Category:Equilibrium chemistry]]
[[Category:Analytical chemistry]]