D.
Brynn Hibbert
^{a} and
Pall
Thordarson
*^{ab}
^{a}School of Chemistry, University of New South Wales, Sydney, NSW 2052, Australia
^{b}The Australian Centre for Nanomedicine and the ARC Centre of Excellence in Convergent Bio-Nano Science and Technology, University of New South Wales, Sydney, NSW 2052, Australia. E-mail: p.thordarson@unsw.edu.au
First published on 25th August 2016
Data analysis is central to understanding phenomena in host–guest chemistry. We describe here recent developments in this field starting with the revelation that the popular Job plot method is inappropriate for most problems in host–guest chemistry and that the focus should instead be on systematically fitting data and testing all reasonable binding models. We then discuss approaches for estimating uncertainties in binding studies using case studies and simulations to highlight key issues. Related to this is the need for ready access to data and transparency in the methodology or software used, and we demonstrate an example a webportal (supramolecular.org) that aims to address this issue. We conclude with a list of best-practice protocols for data analysis in supramolecular chemistry that could easily be translated to other related problems in chemistry including measuring rate constants or drug IC_{50} values.
The three, somewhat linked, key challenges in the analysis of data in supramolecular chemistry concern: (i) determining the stoichiometry of interactions (1:1, 1:2, 2:1,…),^{1} (ii) picking the most appropriate binding model (non-cooperative, cooperative…)^{2} and (iii) obtaining values of thermodynamic quantities such as binding constant(s) K_{i}, with reasonable estimates of measurement uncertainty.^{1–3} It has been stated by one of us that any analysis without proper information on the reliability of results is useless,^{4} which highlights the importance of using well-grounded methods to ensure reliability of the information obtained.
The demand for reliability is now becoming intimately linked to a push for openness and transparency in how the data is handled. Proponents of the Open Science movement have stated its goal as encompassing transparent processes where good practices are characterised by: free, public access to scientific communication, open access to web-based tools that facilitate scientific collaboration and public availability and reusability of data.^{5} This approach has obvious benefits for improving the reliability of results obtained from data analysis. If both the raw data and the methodology used (e.g., software code) is made accessible and as transparent as possible, then it should be easier to detect and correct any mistakes, even post-publication.
This feature article covers some recent developments in supramolecular chemistry data analysis in terms of the three aforementioned challenges with particular focus on uncertainty estimations. The potential role that open science and online tools have in addressing the challenges will also be discussed. First the paper will discuss the related challenges of stoichiometry and selection of a binding model. After brief comments on accuracy in the software used and on method selection (NMR vs. UV-vis), the paper will examine some different approaches to estimating measurement uncertainties using sample data to illustrate these methods. The paper finishes by discussing the role of open science and online tools for data analysis in supramolecular chemistry. To conclude a suggested best-practice list is offered for the researcher. The work here builds on our earlier discussions on the fundamentals of analysing binding data.^{1,2} After a brief review of the key equations for 1:1 and 1:2 equilibria and quality of fit indicators, we shall frequently refer the reader to these papers for more in-depth descriptions of key principles and concepts.
The availability of sufficient computing power to apply more realistic statistical and mathematical approaches has caused a shift from forced assumptions of linearity, pseudo-first-order processes, and graphical methods, to numerical solutions of non-linear systems yielding models and parameters with GUM-compliant uncertainties.^{6} We hope that the current article will assist supramolecular chemists in critically evaluating past and present results and planning their future work to obtain more reliable results.
(1) |
Usually, the concentration of the host–guest complex cannot be obtained directly but it can be related back to the known total concentrations of the host ([H]_{0}) and guest ([G]_{0}) and the equilibrium constant K_{a} through the following quadratic equation:
(2) |
In supramolecular titration experiments, the host is then typically titrated with a solution of guest and the change (ΔY) in some physical quantity that is sensitive to the formation of the host–guest complex is then measured. For UV-vis titration the change in absorbance is measured (i.e. ΔY = ΔA), which is proportional to the concentration of the host–guest complex [HG] from (2) multiplied by the difference between the molar absorptivities of host–guest complex and free host ε_{ΔHG} = ε_{HG} − ε_{H}^{1,2}
ΔY = ΔA = ε_{ΔHG}([HG]) | (3) |
For NMR titrations the observed change ΔY = Δδ is likewise directly proportional to the change δ_{ΔHG} in the NMR resonances between the host–guest complex (δ_{HG}) and free host (δ_{H}) but this time multiplied by the amount fraction of the complex HG:
(4) |
1:2 equilibria can similarly be described through the step-wise equilibrium constants: K_{1} for the formation of 1:1 complex HG, and K_{2} for the formation of 1:2 complex HG:^{1,2}
(5) |
(6) |
Cooperative 1:2 systems are characterised by K_{1} ≠ 4K_{2} whereas for non-cooperative systems, K_{1} = 4K_{2}, simplifying the data analysis as discussed below.
For systems where the measured physical change (ΔY) depends on the amount fraction, such as in NMR titrations, we can use similar approaches as outlined above for simple 1:1 equilibria, to obtain:^{1,2}
(7) |
For UV-vis titrations the right-hand side of (7) is multiplied by [H]_{0}. The concentration of free guest [G] cannot usually be obtained directly but it can be related back to the total concentrations of the host, guest and equilibrium constant through the cubic equation analogue of (2) (not shown here).^{1,2} For other equilibria such as 2:1, similar approaches are then used to obtain relationships between measurable parameters, a physical change that occurs upon forming the host–guest complex, and the equilibrium constants of interest. For 1:1, 1:2 and 2:1 host–guest systems straightforward analytical (exact) solutions are available, but for more complex systems, some shortcuts or approximations are necessary.^{2}
We now turn our attention to the data fitting process. A general model of the problem such as (3) or (7) for N data (x_{data}, y_{data}) with individual observations (x_{i}, y_{i}) can be written as:
y_{i} = f(x_{i},β) + e_{i} | (8) |
(9) |
(10) |
(11) |
(12) |
Picking the correct model is not straightforward. A good fit of a “simple” model does not prove the model as there are usually an (almost) infinite number of other, usually more complicated, models that might fit the experimental data equally well. Traditionally, supramolecular chemists opt for the simplest plausible model (Occam's razor) once the stoichiometry has been determined. However model selection is a mature statistical field and information theory gives approaches (e.g. Bayesian Information Criterion) that could be applied here.^{10} We describe below simple F-value calculations to aid the choice of model.
To narrow down the number of plausible binding models, knowledge about the host–guest stoichiometry is therefore paramount. Once that has been achieved the more subtle differences between available binding models can then be considered.
Recently published work by Jurczak and co-workers^{16} at the Polish Academy of Sciences challenged this view and, in our opinion, essentially spelled the death of the Job plot as a useful tool in analysing supramolecular binding interactions! Their simulations show that the observed maxima in the Job plot (x_{max}) for various cases of 1:2 equilibria may or may not give the “expected” x_{max} ≈ 0.33 for a 1:2 system (Fig. 1). In other words, the observed x_{max} value is often misleading; of all the 1:2 cases shown in columns 3–5 in Fig. 1, only 4 out of 12 have x_{max} ≤ 0.4 and only one is reasonably close x_{max} ≈ 0.33 (0.36 at the bottom right corner of Fig. 1). What Fig. 1 shows is that the Job plot is more sensitive to the K_{1}/K_{2} ratio and host concentration than to the real stoichiometry.
Fig. 1 Simulated Job Plots for various cases of 1:2 stoichiometry with the exception of the 1:1 complexes in column 2 (K_{2} = 0). In all cases K_{1} = 1000 M^{−1} with Y_{HG} = 1 and Y_{HG2} = 2 (additive model). The concentration of host in M is shown in column 1. The simulated maxima in the Job plots is shown as x_{max}. Reproduced with permission from ref. 16. |
Jurczak and co-workers then go further and show that even at reasonably high host concentration (0.01 M) with a low K_{1}/K_{2} ratio of 4, corresponding to classical non-cooperative binding, the outcome is highly dependent on the ratio and direction (increasing/decreasing) of the physical (analytical) signal Y from the 1:1 (Y_{HG}) and 1:2 (Y_{HG2}) complexes formed (Fig. 2).^{16} In NMR titrations, Y_{HG} correspond to the chemical resonance (change) from the 1:1 complex (δ_{HG}) and Y_{HG2} to the one from the 1:2 complex (δ_{HG2}) (see also (7) above).
Fig. 2 Simulated Job Plots for various cases of 1:2 stoichiometry with concentration of host = 0.01 M, K_{1} = 1000 M^{−1} and K_{2} = 250 M^{−1}, i.e. non-cooperative 1:2 binding. Only Y_{HG} and Y_{HG2} between cases A–D. The corresponding case for Y_{HG} = 1 and Y_{HG2} = 2 is shown in row 4, column 5 of Fig. 1. Reproduced with permission from ref. 16. |
These simulations show that for a true 1:2 equilibria, depending on the combinations of Y_{HG} and Y_{HG2}, the observed x_{max} may lie anywhere between 0.29 and 0.63, depending on the assumed 1:2, 1:1 or 2:1 stoichiometry! The message from the data in Fig. 1 and 2 is that Job plots are exceptionally poor indicators of stoichiometry in supramolecular host–guest chemistry. It is for this reason that we propose to declare the Job method as practically dead as an analytical tool in supramolecular chemistry. Jurczak and co-workers point out that the Job method may still have a valid use in the study of inorganic complexes^{16} where the K_{i}'s are typically ≫1/host concentration [H]_{0} (or in other words, the dissociation constant(s) K_{d} ≪ [H]_{0}) which correspond somewhat to the case in the bottom right corner of Fig. 1. This situation is relatively rare in classical supramolecular host–guest chemistry, particularly the type that is studied by NMR titration where the K's hardly exceed 10^{5} M^{−1} for practical reasons.^{1,17} Alternatively, in the case of positive cooperativity, if K_{2} is comparable or even larger than K_{1}, a Job plot will probably give the correct answer (consider this as an extreme case of the ones shown in the right-most column in Fig. 1).
Are Job plots then useful at all? In the limiting cases where either both K_{1} and K_{2} are large or K_{2} is relatively large compared to K_{1}, Job plots still appear to give the “right” answer. However, one can only be certain about this if one has prior knowledge of K_{1} and K_{2}. In other words, at best, a Job plot can only be used for additional positive confirmation of a binding model about which there is sufficient information, i.e. the values of K_{1} and K_{2}. It is practically useless as a tool to rule out more complex models or to differentiate between different binding models. In light of our discussion below, the additional time and effort required to obtain a Job plot, would be much better spent on repeat experiments or by performing a titration experiment at a different concentration of the host.
There are, as outlined in our earlier paper,^{1} other methods than the continuous variation method that can be used to determine stoichiometry in host–guest systems, including the consistency of the model(s) proposed to changes in concentration.^{1} Jurczak and co-workers point out in their paper that a residual plot is probably most useful.^{16} A regular sinusoidal distribution of the residual indicates the assumed model is incorrect but unfortunately, such an observation does not direct the researcher to the correct stoichiometry. In essence, this means that if there any ambiguity about the binding model or the correct stoichiometry, the best approach is to fit the raw data to all probable models and then compare the results.
Fig. 3 (A) The structure of the host (H) 1 and its 1:2 host–guest (HG_{2}) complex 2. (B) The four different binding models (flavours) based on (7) that can be used to describe a 1:2 equilibria. Reproduced with permission from ref. 18. |
The host 1 (Fig. 3A) can bind up to two cations such as Ca^{2+} or Mg^{2+} to form the 1:2 complex 2 (1 can also bind two anions). The binding data for Mg^{2+} and other cations and anions was fitted to all four flavours of the 1:2 equilibria and the results then systematically compared in terms of quality of fit indicators, residual plots vs. number of parameters obtained (Table 1).^{18}
Binding model^{a} | K _{1}/M^{−1} | K _{2}/M^{−1} | cov_{fit} ratio^{b} | SS_{y}/10^{−3}^{c} | df^{d} | F^{e} | P^{f} |
---|---|---|---|---|---|---|---|
a The four different binding models compared (see Fig. 3B). b Raw cov_{fit} from (12) divided by cov_{fit} for the statistical model = 3.66 × 10^{−3}. c Calculated from (10). d Degrees of freedom = N − k. e F-value from (11). In all cases the more complex model (2 in (13)) is the full 1:2 model. f P-value (significance test).^{20} g The SS_{y} for the Noncoop model is 0.004% greater than for the full model. h For noncooperative/statistical binding K_{2} is calculated as K_{2} = K_{1}/4 from the K_{1} value obtained. i The K_{1} value obtained is physically improbable. If as in ref. 18 no constraints are used the model converges on a negative K_{2} which is physically impossible. | |||||||
Full | 4139 | 1059 | 5.6 | 1.98^{g} | 86 | N/A | N/A |
Noncoop. | 4252 | 1063^{h} | 5.6 | 1.98^{g} | 87 | 0.004 | 0.95 |
Additive^{i} | 3 × 10^{6} | 1784 | 0.5 | 7.33 | 90 | 58.1 | 10^{−23} |
Statistical | 15986 | 3997^{h} | 1 | 11.17 | 91 | 79.9 | 10^{−30} |
The most useful indicator used to select a binding model in this study was cov_{fit} obtained from (12). The more complex model, and hence the number of parameters fitted, the better the fit generally is. To justify the selection of a more complex model such as the full 1:2 model over a simpler 1:2 additive model, cov_{fit} needs to be at least three to five times better (lower) for the complex model(s).^{18} For example, in the case of Mg^{2+} binding, the cov_{fit} for both the full and noncooperative 1:2 model was 5.6× better (lower) than the simpler 1:2 statistical model (Table 1). The additive model was ruled out not only because of the poor cov_{fit}, but also as the magnitude of K_{1} obtained was physically improbable. The conclusion from this work was therefore that both the full and noncooperative 1:2 models could be used to describe the binding of Mg^{2+} to 1.^{18}
The above process for selecting binding models relies on the subjective assessment of indicators. A greater number of coefficients will usually achieve a better fit, so a means of taking into account the reduced degrees of freedom would lead to comparable figures of merit. A more statistically robust approach is to test the sum-of-squares from each model by an F test at the appropriate degrees of freedom of each model:^{21}
(13) |
Here, number 1 and 2 refer to the simpler (Noncoop., Additive or Statistical) and more complex (Full) models being compared, respectively, SS_{1} and SS_{2} are the SS_{y} values calculated according to (10) and df_{1} and df_{2} are degrees of freedom calculated from df = N − k with N = number of data points and k = number of parameters. The probability (P) of finding the observed F-value given the null hypothesis that the sums-of-squares are drawn from the same population, i.e. there is no difference between the models, can be readily calculated.^{20} Rejecting the null hypothesis at, say, the 95% level (P < 0.05) implies that the more complex model (number 2 in (13)) does fit the data better than the simpler one (number 1 in eqn (13)). This is not same as saying the more complex model is correct if the P-value is low but that the fit of the data is better described by that model.
We analysed the data shown Table 1 using the F-test and calculated the corresponding P-values (Table 1). The results clearly show that we accept the null hypothesis (P > 0.05) between the full and noncooperative 1:2 binding model, and therefore infer the noncooperative binding model. This is in contrast to the difference between the more complex full model with either the additive or statistical 1:2 model which give minuscule P-values (<10^{−23}). The sum-of-square test yields the same conclusion as the simple semi-subjective quality of fit comparison but it is quantitative and objective.
The frequent use of approximations in older literature and software programs also raises issues. Many legacy programs that are still quite popular use the method of successive approximation to solve the quadratic eqn (4) for 1:1 equilibria or the cubic equation that underpins (7) for 1:2 equilibria. This is no longer necessary as the combination of modern computer processing power and highly sophisticated programs (languages) such as fast and accurate mathematical and statistical algorithms within the open source Python programming language or commercial packages like Matlab solve these polynomials quickly and accurately. Although these legacy programs often get the answer almost right, we have demonstrated previously^{26} that even in the case of relatively simple 1:1 NMR models, the binding constants obtained by one popular legacy program differs by a few percent when compared to a Matlab-based program^{1} that solves (4) directly.
The main limitation of UV-vis spectroscopic titrations is not only the need for a suitable chromophore, but that the titration has to be performed within the absorbance range that follows the Beer–Lambert law, limiting the available concentration range. Fluorescence titration can lower the concentration limit for UV-vis (around 10^{−6} M for strong chromophores), even towards the nM range, but fluorescence titrations are also limited to the relatively narrow concentration range that yields a linear (Beer–Lambert like) response.^{1,2}
Spectroscopic UV-vis or fluorescence titration methods at low concentrations will often “mask” the presence of higher stoichiometries, e.g. in the case of 1:2 host–guest complexes; unless K_{2} is particularly large, the concentration of the 1:2 host–guest (HG_{2}) complex will be minuscule and not detectable in the data obtained.
The “Guide to the expression of uncertainty in measurement” (GUM)^{6} and its supplements is produced by an international collaboration of eight organisations (including IUPAC), the Joint Committee for Guides on Metrology. It offers guidance on how the uncertainty of a measurement result can be estimated from knowledge of systematic and random factors that influence the result. There are strategies to obtain uncertainties that are not provided by statistical treatment of replicate measurements. Taking the measurement equation, uncertainties in each term are combined by the law of propagation of error. Combination of uncertainty terms in quadrature is correct for the linear case, and is also sufficient for problems that are mildly non-linear. However the use of Monte Carlo methods is recommended for obtaining coverage intervals of many problems where non-linearities do not allow simple error propagation.^{29} In analytical chemistry, the elimination, or correction for, systematic errors is a major problem in assuring the quality of results. Differences between results reported by laboratories, often in excess of estimated uncertainties, can be attributed to unknown, and uncontrolled bias.^{30} Error models used to obtain coefficients in supramolecular chemistry always assume the absence of bias, even when no great efforts have been made to demonstrate its absence.
If the 100P% coverage interval is needed (e.g. P = 0.95) it is recommended that considerably more than 1/(1 − P) trials are taken. If M ∼ 10^{4} × 1/(1 − P) the parameters can be expressed to about two significant figures (a relative uncertainty of ∼1%) As each trial to obtain equilibrium parameters requires an iterative, non-linear fit it is not practical with present computing power to run thousands or tens of thousands of trials. It is therefore recommended (Section 7.2.3 of ref. 6)^{6} that if M = 100 or less, and a Gaussian distribution of the parameter values is assumed, then the mean and standard deviation of the set of M parameter values should be used to construct the coverage interval. In the simulations discussed below, M = 200 which means the relative uncertainty on the uncertainty values obtained is closer to 10% than ∼1%, e.g. for a reported uncertainty value below of 8%, the uncertainty on that number is in the order of ±0.8% (rounded to ±1%).
In host–guest titration data analysis one would include the uncertainties of the input concentrations ([H]_{0} and [G]_{0}) and the “best” fitted physical values (y_{calc} in eqn (8)). To estimate the correct variances for these inputs, an uncertainty budget estimation^{3} on the concentration of the solutions prepared and the precision of the observed signal (e.g. chemical shift in NMR) should be made.
(14) |
It is stressed that the standard uncertainty given by (14) has many assumptions, and delivers a symmetrical interval, also called the asymptotic error.^{21}
An alternative approach uses the ‘profile likelihood’ or ‘model comparison’^{21} method in which keeping k − 1 parameters constant at their optimised values, varies the kth to construct a coverage interval of 100 − α% based on a likelihood ratio that is at the α/2 limit of significance, as determined by an F-test using^{21}
(15) |
The value of the binding constant(s) (K) should be the arithmetic average of the N results K_{i}, or, if it is thought that the different data could have significantly different measurement uncertainties, the weighted mean (_{w}).
(16) |
(w_{i} = 1/u_{i}^{2}) | (17) |
The associated combined uncertainty is more difficult to determine. Duewer^{32} gives eleven ways of combining the individual quoted uncertainties and the standard deviation of the means. The two recommended here as practical and easily implemented are the simple sample standard deviation of the mean (s()) calculated from the data with no regard to the individual standard uncertainties (the assumption is the variability between values reflects the internal variability)
(18) |
(19) |
The fit of the data to the full 1:2 binding model was analysed (Table 2 and Fig. 4) using standard uncertainties of the binding constant values (u(K_{i})) from (14) (asymptotic errors), the profile likelihood or model comparison based on (15), and from Monte Carlo simulations (M = 200) based on random sampling from distributions of the input concentration data for the host ([H]_{0}) and guest ([G]_{0}), and the ideal calculated fit data (y_{calc}). The relative standard deviations of the distributions were 2% for [H]_{0}, 1% for [G]_{0} and 0.5% for y_{calc}. The difference in relative standard deviations for [H]_{0} and [G]_{0} can be rationalised based on the fact that the concentration of the guest solution used in supramolecular titration is generally greater than that of the host. Modern NMR (and UV-vis) instruments are also highly accurate making the 0.5% estimation of relative standard deviation conservative.
Binding constant analysed | Type of limit | u(K_{i})^{a} | Model comparison^{b} | Monte Carlo method^{c} | ||
---|---|---|---|---|---|---|
±Limit (%) | Lower limit (%) | Higher limit (%) | Lower limit (%) | Higher limit (%) | ||
a Relative standard uncertainty or asymptotic error^{21} from (14). b Based on (15),^{21} also sometime also referred to as the profile likelihood method (see Fig. 4 for illustration). c From Monte Carlo (M = 200) simulation using 2% relative uncertainty on [H]_{0}, 1% relative uncertainty on [G]_{0} and 0.5% relative uncertainty on y_{calc}. The uncertainty values obtained from Monte Carlo have themselves approximately ±10% relative uncertainty. d Standard deviation = standard uncertainty according to (16) (calculated using de Levie's method)^{31} or 68.3% confidence interval (P-value or α = 0.317) for the Model Comparison and Monte Carlo methods. e The 95% coverage interval (CI with P-value or α = 0.05). For the standard uncertainty method, the value obtained from (16) is multiplied by the t-value at α = 0.05 and divided by √N. | ||||||
K _{1} | u ^{ } | ±4.1 | −12 | 14 | −9 | 9 |
U _{95%} ^{ } | ±8.2 | −15 | 18 | −16 | 24 | |
K _{2} | u ^{ } | ±5.2 | −14 | 18 | −4 | 3 |
U _{95%} ^{ } | ±10 | −19 | 24 | −7 | 6 |
Fig. 4 Graphical representation of the data shown in Table 2. The diamonds represent one of the M = 200 simulated results from the Monte Carlo calculations. The 68.3% or one standard deviation (blue) and 95% (black) confidence intervals for both the Monte Carlo (broken thick lines) and Model Comparison (solid thin lines) or Profile Likelihood methods. The symmetrical standard uncertainty (inner dotted red line box) or asymptotic error from (15) (calculated using de Levie's method)^{31} limits are also shown together with the corresponding 95% confidence interval on the standard uncertainty value (outer dotted red line box). The insert shows two of the outliers (arrows) from the Monte Carlo simulations. |
The differences between these methods are most readily seen when they are plotted graphically as a function of relative (%) error of K_{1} and K_{2} (Fig. 4). The scatter pattern that the 200 Monte Carlo simulations, and the corresponding confidence limits (broken black and blue lines) are clearly not symmetrical around the “best fit” K_{1} and K_{2} values (0%). The Monte Carlo results show more spread or up to 50% from the best fit K_{1} value along the K_{1} axis, but at the most only about 10% away from best fit K_{2} value. The Model Comparison (Profile Likelihood) method is not symmetrical either but has a distinctly different shape than the Monte Carlo with the uncertainty being larger along the K_{2} axis. In contrast, the standard uncertainty or asymmetric error is symmetrical, unlike the raw Monte Carlo scatter.
These results demonstrate the approximate nature of estimating standard uncertainty by eqn (14) or asymptotic error methods. Monte Carlo and Model Comparison (Profile likelihood) methods also give quite different results. The scatter in the Monte Carlo simulations suggests there is better information (smaller uncertainty) on K_{2} than K_{1} in this system. This makes good sense; the ratio of K_{1}/K_{2} suggests noncooperative binding (see also Table 1) and the calculated amount fractions (see unique URL from supramolecular.org^{19} and Fig. S3 in ESI†) shows that the maximum amount fraction for the formation of the 1:1 complex is 0.5, where the 1:2 complex reaches a amount fraction of 0.9 at the end of the titration. This means the NMR data obtained from this experiment “sees” the 1:2 complex better than the 1:1 complex, resulting in a smaller uncertainty on K_{2}. In line with best practice recommendation in the GUM,^{29} the Monte Carlo results appear to give the best presentation of the underlying uncertainties in host–guest binding studies.
Encouraged by these results, we carried out a large number of Monte Carlo simulations for NMR titrations for both the 1:1 and 1:2 binding equilibria, using a range of binding constants (K_{i}). In all cases the chemical shifts are assumed to be additive (see Fig. 3B), the total host concentration was fixed at [H]_{0} = 10^{−3} M which is typical for NMR titrations and the final guest concentration at [G]_{0} = 0.035 M (35 equivalents). The results for the 1:1 binding equilibria are shown in Fig. 5.
The results give valuable insight into factors that affect uncertainty in determining binding constants from host–guest titrations with some interesting but expected trends clearly evident. In terms of the lower limits (Fig. 5, left panels) and the higher limits (Fig. 5, right panel) on the uncertainties on the expected K_{a}, they are not quite symmetrical with slightly larger uncertainties on the higher limits for a given relative standard deviation of [H]_{0}, [G]_{0} and y_{calc} (δ_{calc}). For K_{a} > 5 × 10^{3} M^{−1} the uncertainties do not exceed 10% regardless of the variance. As the K_{a} get larger, the variance has a more pronounced effect on the estimated uncertainty. At K_{a} > 10^{5} M^{−1}, which in any case is close to the limit achievable by NMR,^{1,17} meaningful K_{a} (uncertainty <40%) estimate can only be obtained with the noise (errors) on [H]_{0}, [G]_{0} and y_{calc} (δ_{calc}) are vanishingly small or 0%.
The results from our simulations for NMR titrations for 1:2 equilibria are shown in Fig. 6. These simulations were carried out under conditions where K_{2} = 0.1 × K_{1} or in other words, mild negative cooperativity (interaction parameter^{12,33}α = 0.4). These simulations start by assuming that the chemical shift changes are additive (δ_{HG} = 0.5 ppm and δ_{HG2} = 1 ppm).
Fig. 6 Monte Carlo simulations (M = 200) for NMR titrations with underlying 1:2 equilibria (full model). In all cases [H]_{0} = 10^{−3} M, [G]_{0} is spread unevenly across 49 data points between [G]_{0} = 0–0.035 M, with δ_{HG} = 0.5 ppm and δ_{HG2} = 1 ppm (as in the additive model – Fig. 3B) as the starting values in the “ideal” datasets used at the start of these simulations. The data was fitted with K_{1} between 10^{2}–10^{6} M^{−1} and K_{2} between 10–10^{5} M^{−1} with K_{2} always fixed at K_{2} = 0.1 × K_{1} (mild negative cooperativity). The variance is either 0% or 1% on y_{calc} (δ_{calc}), 0–2% on [G]_{0} and 0–4% on [H]_{0}. In each case the variance of [H]_{0} is 2× that of [G]_{0}. The contour plot is coloured according to relative (%) the Monte Carlo uncertainty on the expected K_{1} (left column) and K_{2} (right column) values at the 95% confidence interval level (95% CI). (A) The colour scheme used in the contour plots between −100% and +100% of the expected K_{1} or K_{2} values. (B) Contour plot of the calculated Monte Carlo 95% CI on K_{1} assuming uncertainty of 0% for y_{calc} (δ_{calc}). The lines indicate steps of 10%. The ±20% (blue) ±40% (orange) and ±60% (red) levels are highlighted with bold lines and labels. The y-axis represents changes in K_{1}. Starting from the central vertical line, the x-axis represents increasing noises in [H]_{0} and [G]_{0} for the lower limit (left side panel) and higher limit (right side panel) of the expected K_{1}. (C) Same as (B) except for K_{2} instead of K_{1} including the y-axis which shows K_{2}. (D) Same as (B) except with additional 1% relative standard deviation on y_{calc} (δ_{calc}). (E) Same as (D) except for K_{2} instead of K_{1} including the y-axis which shows K_{2}. Note the order of magnitude difference in the y-axis between the left (K_{1}) and right (K_{2}) columns. |
These results immediately demonstrate that the estimated uncertainty is highly sensitive to the relative standard deviation on both the total host and guest concentrations and the calculated y_{calc} (δ_{calc}) values. The higher limit on K_{1} is particularly vulnerable to any variance; if the relative standard deviation of [G]_{0} exceeds 0.5 (relative standard deviation on [H]_{0} < 1%), a meaningful estimate on the upper limit of K_{1} cannot be obtained (Fig. 6B and D). The uncertainty of K_{1} is on the other hand fairly insensitive to variance on the calculated y_{calc} (δ_{calc}) values, or the expected measured analytical signal. This suggests that researchers need to take particular care at preparing host–guest solutions when a 1:2 equilibria is suspected to minimize the resulting uncertainty.
Interestingly, the simulations also suggest a region of stability where one expects fairly accurate K_{2} values, i.e. when K_{2} > 10^{3} M^{−1}. In contrast if K_{2} < 10^{3} M^{−1}, it appears that it would be very hard to obtain a meaningful estimate on K_{2}. This does make sense as with a low K_{2} there would little “information” about the 1:2 complex on the expected binding isotherm whereas for high K_{2} the opposite is true (a high K_{2} and a 10× higher K_{1} would be beyond the practical limit for NMR).^{1,17}
Up to this point we have only discussed how the uncertainty on binding constant(s) obtained from host–guest titration studies could be estimated from a single experiment (n = 1). This is, however, not the ideal situation – if possible, an experiment needs to be repeated several times and the uncertainty then estimated from these n-repeats. The uncertainty of each individual fit then becomes incorporated in the uncertainty or variations across the multiple repeat, allowing us, as we show below, to largely ignore the estimated uncertainty on individual fits. At the very minimum, one should perform the experiment in triplicate (n = 3) but we strongly recommend that titration experiments should be carried out in quadruplicate (n = 4) as the results for n = 4 are significantly more statistically robust than for n = 3. Compared with three repeats, four repeats will improve the ratio of the t-values multiplied with the ratio of the square root of n (see (18)) by 36% and going from n = 4 to n = 5 improves that number by 21%.
The next question is how to estimate the uncertainty on the arithmetic or weighted mean values we obtain. To answer this we return to our Mg^{2+} titration of 1 example (Fig. 3 and Table 1) but this time looking at three repeats of this experiments (n = 3) with the results summarised in Table 3.
K _{1} value/M^{−1} | Uncertainty^{a} | K _{2} value/M^{−1} | Uncertainty^{a} | |||
---|---|---|---|---|---|---|
−s^{b}/−M^{−1} | +s^{c}/M^{−1} | −s^{b}/M^{−1} | +s^{c}/M^{−1} | |||
a Standard deviation (68.3% confidence interval, CI) limit on K_{i} obtained from Monte Carlo (parameters the same as in Table 2 and as M = 200, these uncertainty values themselves also have a relatively uncertainty of ±10%). b The lower 68.3% CI limit from Monte Carlo. c The upper 68.3% CI limit from Monte Carlo. d Mean = sum of K_{i} values/n. e Calculated standard deviation of the mean from (18). f Weighted mean from (16); lower value calculated from lower Monte Carlo limits (−s) and higher value from the higher Monte Carlo limits (+s). g Weighted standard deviation from the mean calculated with (19) and based on the lower (−s) and higher (+s) Monte Carlo uncertainties. | ||||||
Experiment 1 | 4139 | −378 | 386 | 1059 | −40 | 30 |
Experiment 2 | 4832 | −480 | 545 | 907 | −38 | 32 |
Experiment 3 | 3840 | −372 | 490 | 667 | −28 | 32 |
(mean)^{d} | 4271 M^{−1} | 878 M^{−1} | ||||
±s()^{e} | ±293 M^{−1} | ±114 M^{−1} | ||||
(Relative s(), %) | (±6.9%) | (±13%) | ||||
_{w} ^{ } | 4185 M^{−1}–4216 M^{−1} | 823 M^{−1}–1050 M^{−1} | ||||
±s(_{w})^{g} | −275 M^{−1}, +261 M^{−1} | −122 M^{−1}, +114 M^{−1} | ||||
(Relative s(_{w}), %) | (−6.6%, +6.2%) | (−14.9%, +10.9%) |
Apart from calculating the (normal) mean and (normal) standard deviation of the mean s() from (19), we also include here the weighted mean _{w}(16) and weighted standard deviation of the mean s(_{w}) (19). And as the Monte Carlo uncertainties are not symmetrical, this result in different lower and upper limit estimated of _{w} and s(_{w}). Interestingly, the weighted _{w} and s(_{w}) do not seem to differ much from the normal and s() values even though only the former incorporates the uncertainty (here Monte Carlo) estimates from the individual fits. We conclude from this that in most cases, the normal (unweighted) process for calculating the mean and standard deviation of the mean is perfectly acceptable.
There are two key lessons to take from this last example; firstly, by doing n-repeats and assigning the (normal) standard deviation of the mean of these repeats to the standard uncertainty, one can almost ignore the problem of estimating the uncertainty on the parameters in each experiment! We put the qualifier almost here, as it would still be prudent to estimate these, even if it just done with the simple standard uncertainty (asymptotic error) method using (15) to check how realistic the obtained K_{i}'s are. The other interesting lesson is that for this experiment, the relative standard deviation of the mean is greater (13%) for K_{2} than for K_{1} (6.9%) – opposite from the estimation of the uncertainties from the individual fits according to Tables 2 and 3! The reason for this becomes evident if one compares the results from Experiment 3 to the others in Table 3 as the K_{2} value appears to be an outlier resulting in a large standard deviation of the mean for K_{2}.
Open access to publications is now quite common but researchers in host–guest chemistry, and to some extent in chemistry at large, have been relatively slow to adopt other important Open Science tools. Depositing raw data with manuscripts is the exception and until now, there has not been any open access database for host–guest complexation data. Contrast this with single crystal X-ray crystallography where for 50 years the Cambridge Crystallographic Data Centre (CCDC) has allowed researchers both to deposit raw data and then search and retrieve any other deposited data for further analysis.^{34}
One of us (Thordarson) has now established the web-portal OpenDataFit^{35} which includes the site supramolecular.org.^{19} This site provides data deposition and storage, and offers the community a free (open) access web-tool to fit their data to a range of binding models. The software code is open source (Python) and available online for scrutiny and further improvements (Fig. S1, ESI†).
The website is built around the concept of end-users uploading their input data (host and guest concentrations and the measured physical signal such as NMR resonance) in a simple spreadsheet format. The user then selects between various binding models and sets parameters such as initial guesses of the binding constant(s) being sought. After fitting the data and examining the results, which include residual and mole-fraction (speciation) plots, the user either refines the binding process further or saves the results. The archiving step includes an opportunity to add metadata such as which host and guest were used, solvent, temperature and other useful information. The user is then given a unique URL that can be used to access the data later as well as the option of downloading all the results in a spreadsheet for further analysis and plotting (Fig. 7 and Fig. S1–S4, ESI†).
Fig. 7 The Open Access website supramolecular.org^{19} for data analysis and archiving in host–guest supramolecular chemistry. (A) Flow diagram of how data is processed on the website. (B) Snapshot from the website showing the result window from data archived at the unique url: http://app.supramolecular.org/bindfit/view/8a658114-0b28-4c63-92c0-09a7a976f0be which was fitted to 1:2 NMR binding data. Additional screenshots are provided in ESI† (Fig. S2–S5). |
The website offers inter alia, global^{36} fitting of host–guest titration data to 1:1, 1:2 and 2:1 binding models, including the various 1:2 (and corresponding 2:1) flavours mentioned above. Users can also choose between the robust Nelder–Mead^{37} (Simplex) algorithm and the L-BFGS-B^{38,39} (limited (L) memory quasi-Newton Broyden–Fletcher–Goldfarb–Shanno simple box (B) constraints), which allows constraint of the search space. The site also returns the standard uncertainty^{21} (asymptotic error – see above) of the estimated binding constants. Users can also simulate binding data to help design their experiments. Planned additions include Monte Carlo estimation of uncertainty (see above). For the 1:1, 1:2 and 2:1 binding models, the fitting processes are based on the exact solutions (no approximations) of the equations used to describe these binding equilibria, e.g., (4) in the case of the 1:1 NMR model.
The database function should in time, allow for systematic investigation of inter- and intra-laboratory biases by comparing data from different laboratories. Data mining might also allow investigators to do more systematic investigations into methods for comparing models, answering questions such as how robust the F-value (13) calculation method really is with a large set of real data.
It is our hope that the supramolecular.org website might also provide a catalyst for the supramolecular chemistry community for setting minimum standards for publishing data, along the lines of what is already common practice in the crystallography community. Other data-intensive fields have also established community-based “minimum information criteria” for publishing data, with proteomics being a prime example.^{40}
(2) If more than one stoichiometry or binding model is suspected, fit the data to all models and systematically compare results to eliminate those that do not fit based on criteria's such as shape (scatter) of the residual plot^{16} and the sum-of-squares F-test^{19} (eqn (13)).
(3) If certain binding constants in multi-species equilibria are very low, e.g. K_{2} in the 1:2 host–guest equilibria, the information content associated with that complex is inherently very limited. No method, no matter how sophisticated is likely to yield any reliable estimate of that binding constant (they will have large uncertainties).
(4) Reviewers should without hesitation request for a revision of any papers that still use outdated and inaccurate linear-transformation methods such as Benesi–Hildebrand.
(5) When possible, the program used to fit the data should calculate the concentrations of the species involved using the exact mathematical expression for the equilibria of interest. For simple 1:1, 1:2 and 2:1 binding model this should be mandatory and legacy programs that use the method of successive approximation should not be allowed.
(6) Best practice for estimating uncertainties clearly involves repeating the experiment at least 3, ideally 4 or more, times. The estimated uncertainty on the fit in individual experiments should be checked for signs of very poor fit. This can be done with the simple standard uncertainty method (asymptotic error) using (14) but Monte Carlo is desirable.^{42} If the individual fits are not unreasonable, they can subsequently be ignored and the uncertainties of binding constants and other parameter then simply estimated from the standard deviation of the mean s() from the n-repeat experiment.
(7) If n-repeat measurements cannot be performed for practical reasons, we strongly recommend that the estimation of the uncertainty on the fitted parameters should, as recommended by the GUM Supplement 1,^{6,29} be performed by Monte Carlo simulations^{42} and reported at the 95% confidence interval level.
(8) Data should be made accessible (Open Access)^{5} and the software code used available and transparent to ensure others can verify and analyse further the experimental data.^{19}
Repeating once again the message that “data without any information about its reliability is meaningless”^{4} we hope that this article and the list above will provide the host–guest supramolecular chemistry community with the necessary tools to make their results even more reliable. This will reduce mistakes and confusion in the field and accelerate further the rapid progress of host–guest supramolecular chemistry. We note also that the above approach could easily be adopted in other fields of chemistry, e.g. when determining rate constants or measuring drug potency via experiments to measure IC_{50} or EC_{50}.
Footnote |
† Electronic supplementary information (ESI) available: Unique URL's for key data, Matlab m-files for Uncertainty Estimations, additional website screenshots. See DOI: 10.1039/c6cc03888c |
This journal is © The Royal Society of Chemistry 2016 |