Analysing binding stoichiometries in NMR titration experiments using Monte Carlo simulation and resampling techniques

Monte Carlo simulation (MC)

The usage of Monte Carlo simulation (MC) and Bootstrapping (BS) to estimate confidence intervals using the percentile methods as proposed previously (Thordarson, 2011; Lowe, Pfeffer & Thordarson, 2012) has already been presented in an earlier article (Hübler, 2022b). Recently, Kazmierczak, Chew & Vander Griend (2022b) have analysed BS in context of photometric titration experiments, including stock solution errors.

After performing MC or BS, the confidence interval does not contain information of the distribution of individual obtained values for the parameters as only the percentile is calculated (Efron, 1979). The dimension of the distribution e.g., the number of Monte Carlo runs and the individual probability of any value to be obtained after a single run are not taken into account. To include the number of Monte Carlo steps and the probability of an individual parameter value, the histograms are further characterised using the standard deviation (Eq. (5)) and the Shannon entropy (Eq. (6)). (5)${\displaystyle {\sigma }_{dist}=\sqrt{\frac{\sum _i^N{\left(x_i-\overline{x}\right)}^2}{N-1}}}$ (6)${\displaystyle H\left(x\right)=-\int p\left(x_i\right)ld\left(p\left(x_i\right)\right)}$

While the standard deviation is a common tool in every-days statistical judgement, the Shannon entropy is more often used in pattern recognition and information theory (Bishop & Nasrabadi, 2006). It is out of the scope of this article and SupraFit to deal with the various aspects in pattern recognition, but the most important facts shall be summarised. If only one possible result for a parameter is obtained, the probability p(x) is 1, resulting in an entropy of 0. If however the possible results for a parameter are uniformly distributed across N possible values, the entropy is maximum with −ln(1/N). As consequence, the lower values of the Shannon entropy indicate sharper peaks, while higher entropy indicate broader peaks. In case of a normal distribution, the standard deviation and the Shannon entropy are directly linked by Eq. (7) and the Shannon entropy gives no additional information. However, for non-normally distributed random numbers (Bishop & Nasrabadi, 2006) like in case of the distribution of the parameter values after Monte Carlo simulation in connection with nonlinear regression, Eq. (7) is not true and the Shannon entropy accounts for the distribution of the random numbers. (7)${\displaystyle H\left(x\right)=\frac{1}{2}\left(1+ln\left(2\pi {\sigma }^2\right)\right)}$

As the Shannon entropy does not depend on the concrete values of the parameters (x) itself but on the probability (px(x)), the Shannon entropy is dimensionless. The standard deviation however has the dimension of the value. The comparison of the results of Monte Carlo simulation with parameters of different dimensions is hence only meaningful if the dimensionless Shannon entropy is used.

The discrete calculation of the Shannon entropy is given in Eq. (8), where B denotes the number of bins and Δ the width of a single bin. (8)${\displaystyle H\left(x\right)=-\sum _i^{B}p\left(x_i\right)\Delta ld\left(p\left(x_i\right)\right)-ld\Delta }$

Since the Shannon entropy is in SupraFit calculated with the same number of bins for each distribution individually, the width of a single bin depends on the distribution of each parameter itself. For more compactly distributed values of a parameter, smaller bins are obtained resulting in artificially increased entropy values. Hence the second term in equation is omitted in the standard calculation and the Shannon entropy is calculated according to Eq. (9) if not requested otherwise.² (9)${\displaystyle H\left(x\right)=-\sum _i^{B}p\left(x_i\right)\Delta ld\left(p\left(x_i\right)\right)}$

Cross Validation (CV)

Cross Validation defines a resampling method (Efron, 1979), which usually is applied to determine the ideal parameters or latent variables in multivariate statistics such as PCA or PLS for example in QSAR studies (Gramatica, 2007). The parameters or latent variables are judged according their predictive behaviour, which can in the simplest case be determined as follows: Having a data set with N elements, in Leave-One-Out Cross Validation (L1O-CV) N-1 data points are used to train a new model using a set of chosen parameters or number of latent variables. The performance of the newly defined model can be evaluated using the omitted data point which forms the test set. As there are N independent possibilities to form a training set, N independent evaluations can be used to judge the performance of the variables. Analogously, in Leave-Two-Out Cross Validation (L2O-CV) pairs of two points configure the test set. The training set is composed of the remaining N-2 data points, with N*(N-1) possible combinations. As the order is not relevant, only N*(N-1)/2 combinations are unique.

In SupraFit, CV is not used in the described manner, which is to recalculate the test set data points and estimate the performance of the variables on the base of statistical judgment. Instead, it is implemented to judge the models if hypothetically fewer data points were acquired. In the simplest approach (L1O-CV), each data point of the experimental data set is left out once and the model is fitted. In the limit of only one tracked NMR signal or an ITC experiment, only single points can be left out. Having more than one NMR signal analysed or several wavelength in case of photometric experiments, each point on the x axis is assigned a vector of observed values, which will be simultaneously left out. However, they will still be referred as data points within this context. As there are again N ways to leave one point out, a distribution of N model parameters can be obtained. The analysis of the distribution is performed in the same fashion as introduced in case of the Monte Carlo simulation. Furthermore, L2O-CV as well as the generalisation Leave-X-Out Cross Validation (with X < N) can be performed and analysed within SupraFit. In case of the Leave-X-Out Cross Validation, the number of trials can be calculated using the general formula for combinations (Eq. (10)). (10)${\displaystyle S_{max}=\frac{N!}{X!\left(N-X\right)!}}$

As advantage over the Monte Carlo Simulation, where an additional input parameter is required and the results are not exactly reproducible, Cross Validation works without the additional input and leads to reproducible results for the same combinations.

The individual parameter fitting procedures during Cross Validation are performed as if there were Monte Carlo simulation, including the parallelisation. A benchmark on the MC implementation with respect to the parallelisation was given previously (Hübler, 2022b). However some aspects of the data preparation, that is (a) the generation of the map containing the pseudo-experimental data with left out points and (b) the selection of the pseudo-experimental data, will be discussed in this context. In general, the valid combination of points to be left out will be precalculated in order to maintain unique calculation during parallel execution.³ In case of L1O and L2O, the number of parameter fitting processes are N and N*(N-1)/2. Both can easily generated using a simple loop or two nested loops, and if N < 100 routinely be evaluated in SupraFit. In case of more points X to be left out and/or significantly larger amount of data points N acquired, the number of possible combinations S_max easily reaches the limit of computational time and resources.⁴ As consequence a random subset of the combinations has to be chosen prior to the parameter fitting. In case of a greater number of N and the chosen value of X, the preparation of all combinations S_max may take longer than the fitting of the subset with S combinations itself. Therefore two approaches for the generation of the map are implemented in SupraFit. (a) The whole map is calculated using nested loops. In case there are fewer steps S than possible steps S_max requested, the individual pseudo-experimental data sets are randomly chosen from the fully generated map using uniformly distributed random numbers. (b) In the second approach, uniformly distributed random numbers are used to generate S sequences (having the length X) of unique numbers which define the points to be left out. Uniqueness of the sequences itself is ensured after sorting a newly obtained one first and then comparing it to the previously generated sequences. A new unique number within the sequences becomes more improbable as the sequence gets more complete (1/X for the last number) and the uniqueness of the sequences becomes more improbable as the number of already stored sequences increases. Therefore, this approach becomes only more efficient compared to the nested loop precalculation if significantly fewer steps S than S_max are requested. Depending on the number S and the ratio of S and S_max, SupraFit choses one approach in favour of the other, but each of the two algorithm can be enforced from the user interface.

Reduction Analysis (RA)

During experimental design, the number of data points to be acquired and the ratio of the initial concentration of the substances B and A to be analysed is of particular interest. It has already been discussed by Thordarson, that for many systems ten points between a ratio $\frac{B_0}{A_0}$ of 0 and 1.5 are most important, and additionally ten more up to a ratio of 50 should be acquired (Thordarson, 2011). SupraFit provides a new protocol called Reduction Analysis (RA) which exploits possible redundancy if data points are included in the data set, which itself do not change the model parameters. The fitted parameter of a 1:1 model on a hypothetical experiment with 1:1 stoichiometry would not differ if fewer data points where included during the fitting process, as long as more data points as essentially necessary for an adequate determination for the 1:1 species are available. This approach is comparable to the Mole-Ratio method: After the saturation point($\frac{B_0}{A_0}>1$) is reached, the observed chemical shift changes less than in the range of ($0<\frac{B_0}{A_0}<1$). Limitation itself may arises if the stability constants are low and saturation of the complex formation is not reached within that threshold ratio, that is defined by the stoichiometry of the complexes. While the condition of sufficiently stable complexes can not be expected a priory, the potential of Reduction Analysis to indicate if model parameters are not appropriate will be discussed within this article.

The automatic procedure is based on the step-wise reduction of the number of data points beginning at the end. After each removal, the model is fitted to the remaining points of the data set. The best-fit parameters $\hat{\theta }$ are stored and finally plotted as “function” of the highest ratio, which is the ratio of the last available data point. Since titration should reach at least the saturation point, all fitted parameter below the particular ratio are assumed to differ significantly, even if the stoichiometric model is correct. To allow a rational comparison of Reduction Analysis performed on different models, the standard deviation of the parameters are calculated. Since every model needs a different ratio to be reached during titration, the comparison of the naive standard deviation is not meaningful. Instead the partial standard deviation σ_pt (Eq. (11)) is calculated, taken all data points above a cut-off ratio into account. This cut-off is defined by the “highest necessary” saturation point, therefore ensuring, that for all models the saturation point has been reached. In SupraFit this cut-off is automatically set for NMR titration experiments according to the tested models but may be changed to any value if necessary. (11)${\displaystyle {\sigma }_{pt}=\sqrt{\frac{\sum _{i=k}^N{\left({\hat{\theta }}_i-\overline{\hat{\theta }}\right)}^2}{N-k-1}}}$ ${\displaystyle \text{with:}\text{N - last data point}}$ (10)${\displaystyle \text{k - first data point}}$

Data set with 1:1/1:2 stoichiometry

Monte Carlo Simulation

Exemplary results of the Monte Carlo simulation (with S = 2000, σ_mc = SE_y) of all four models on top of the 1:1/1:2 data set are given in Figs. 1 and 2. The plots show both, all parameters and stability constants only. The factor analysis of the chemical shifts using Sivvu.org reveals three factors (Fig. S20B), indicating that three species contribute to observed signals. This is in agreement with the original model, however no further information about the stoichiometry of the species can be deduced using the analysis.

Following the “visual inspection,” the comparison of the histograms already gives an estimation of the quality of the used models. Since the 1:1/1:2 model recovers the original parameters, the results of the Monte Carlo simulations on top of the best-fit parameters of the corresponding 1:1/1:2 model are considered as ideal results, having appropriate distributions of the parameters and therefore standard deviation and Shannon entropy values indicating well estimated model parameters. The histogram of each parameter shows a narrow distribution (Fig. 1), which is also true for the simpler 1:1 model. The MC simulation results of the 2:1/1:1 model give rather broad distributions (Fig. 2a), and since all parameter are plotted at once, they overlap. Neglecting chemical shifts and visualising only the stability constants, the histograms reveal a broader distribution of the lg K₂₁ and a narrower distribution of the lg K₁₁ values (Fig. 2B). The values of the stability constant remain physical meaningful between one and three. The histogram of the most complex model show only broadly distributed parameters, and focusing on the stability constants only, lg K₂₁ ranges from −7.5 to 3.0 being clearly non-physical. The correct parameter appear to be narrowly distributed. Hence, from pure “visual inspection” the 1:1/1:2 model and 1:1 model with both having narrower distributed parameters appear to behave much better than the 2:1/1:1 and 2:1/1:1/1:2 model with broader distributed parameters.

However, the “visual inspection” is rather biased and need strictly comparable charts (scaling, resolution, etc.). A more unbiased way utilises statistical descriptors derived from the distribution of the parameters itself. In case of the stability constants, these statistical descriptors—the confidence intervals and the values for the Shannon entropy and standard deviation—are given for all four tested models in Table 1. As stated above, the results from the 1:1/1:2 model are taken as reference. The results of the remaining models will then be compared to the reference result.

Table 1:

Stability constants of the four tested models, the 95% confidence intervals obtained from Monte Carlo simulations including the median, σ_dist and H(x) calculated from the distribution of the parameters, partitioned into 30 bins.

parameter	$\hat{\theta }$	${}_{\Delta \theta -}^{\Delta \theta +}$	[θ−	θ+]	σdist	H(x)
1:1/1:2 model; SSE = 0.0001; SEy= 0.0011
lgK11	3.81	+0.04 −0.04	3.77	3.85	0.0188	0.0523
lgK12	2.11	+0.04 −0.05	2.06	2.15	0.0242	0.0640
1:1 model; SSE = 0.0365; SEy= 0.0171
lgK11	3.10	+0.07 −0.07	3.03	3.17	0.0344	0.0837
2:1/1:1/1:2 model; SSE = 0.0001; SEy= 0.0011
lgK21	1.99	+0.51 −4.53	−2.54	2.50	1.6933	1.4813
lgK11	3.81	+0.06 −0.06	3.75	3.87	0.0294	0.0796
lgK12	2.04	+0.09 −0.10	1.94	2.14	0.0496	0.1164
2:1/1:1 model; SSE = 0.0018; SEy= 0.0039
lgK21	1.69	+0.49 −0.50	1.18	2.18	0.2637	0.4507
lgK11	2.67	+0.09 −0.12	2.55	2.76	0.0505	0.1336

DOI: 10.7717/peerjachem.23/table-1

All parameters are unsymmetrically distributed around the best-fit values, which is in agreement with the general findings by Motulsky and Christopoulos (Motulsky & Christopoulos, 2003). A detailed analysis in case of stability constants has later been done by Thordarson (Thordarson, 2011) and Vander Griend (Kazmierczak, Chew & Vander Griend, 2022b). In the 1:1 model, the confidence limits for lg K₁₁ are a nearly twice the limits of lg K₁₁ in the 1:1/1:2 model. Furthermore, the standard deviation and Shannon entropy for lg K₁₁ are higher in the pure 1:1 model. According to the descriptors, the correct parameters in the 2:1/1:1/1:2 model behave slightly worse than in the 1:1/1:2 model, but the confidence interval as well as σ_dist and H(x) are in the same order of magnitude. On the other hand, for the parameter lg K₂₁ the simulation lead to the widest confidence interval and the worst values of σ_dist and H(x). The parameter lg K₁₁ in the 2:1/1:1 model is connected to a value of the Shannon entropy between one obtained in case of the pure 1:1 and the mixed 2:1/1:1/1:2 model, but with simultaneously having the worst values of σ_dist and largest confidence interval. In contrast to lg K₂₁ in the 2:1/1:1/1:2 model, the statistical outcome in the 2:1/1:1 model is much better: a tighter confidence interval and smaller values of σ_dist and H(x) were obtained. As lg K₁₁ and lg K₁₂ are sufficient to recover the simulated titration curve, lg K₂₁ can freely accept any value with only slightly improving the resulting SSE. If the parameters linked to the chemical shift were not included in charts, from pure inspection of the stability constant without comparing to other models, the 2:1/1:1 model could still be acceptable, with the unknown source of broader confidence interval maybe rooted in experimental circumstances. Taking all parameters into account (Figs. 2A and 2C), it is without comparing to other models clear, that the 2:1/1:1 and the 2:1/1:1/1:2 model are unsuited for the data set since some of the individual parameter take even negative values. On the other hand using this judgment, the 1:1 model describes the data sufficiently well.

Figure 3 shows the Shannon entropy (H(x)) and the standard deviation (σ_dist) of the Monte Carlo simulations with 2000 steps compared to simulations with 10000 steps. The qualitative behaviour of both descriptors is identical for the tested models: The lowest obtained values are assigned to the parameters of the 1:1/1:2 model. According to the individual values of descriptors characterising the stability constants, lg K₂₁ is inappropriate while lg K₁₁ and lg K₁₂ are adequate to obtained well fitted models.

The former analysis indicates, that in order to ensure the qualitative behaviour of the statistical parameters, the number of steps S in Monte Carlo simulations is of less importance. However, for each model the according σ_MC value obtained after fitting was used, which is the default input for the Monte Carlo simulation. Hence the different outcome of the Monte Carlo simulation could be traced back to the different σ_MC values. In a subsequent analysis, the input standard deviation σ_MC was then varied as follows: Monte Carlo simulations of each best-fit model were repeated with all possible SE_y value obtained from the best-fit parameters of the four models. The results, visualised as bar charts in Fig. 4, show the dependency of the Shannon entropy from the input standard deviation for each model. In Fig. 4, results depicted with the letter A were obtained using SE_y from the 1:1 model (SE_y = 0.01708) and with the letter B SE_y from the 2:1/1:1 model (SE_y = 0.003889) was used. The results depicted with C are linked to the 1:1/1:2 model (SE_y = 0.00106) and with D to the 2:1/1:1/1:2 model (SE_y = 0.00108).

The 1:1 model with the fewest parameters exhibits the worst fit, which is indicated by the highest value of SSE and SE_y (compare Table 1). Hence the individual values of the descriptors obtained from Monte Carlo simulation with input A indicate worse performance than in case of the input B to D, which are connected to lower values of SE_y. The bar charts in Fig. 4 show clearly lower values of H(x) in case the input SE_y was taken from a more complex model. Furthermore, the entropy falls below the values obtained using correct 1:1/1:2 model when the lower SE_y from the 2:1/1:1/1:2 model is applied. The SE_y value from the 1:1/1:2 model (C) is the second-best result, therefore using the SE_y value obtained from the fit of the 1:1 model (A) or 2:1/1:1 model (B) worsen the results of the Monte Carlo simulations clearly. As the SE_y obtained from the 2:1/1:1/1:2 model (D) is the lowest, the best σ_dist and H(x) results are obtained using this input. In the same manner, the resulting entropy gets lower when the Monte Carlo simulations of the 2:1/1:1 model are performed with SE_y input from 1:1/1:2 or 2:1/1:1/1:2 model (Fig. 4B). On the other hand, the Shannon entropy increases if the SE_y value from the 1:1 model is applied. Common for all calculations is the bad performance of lg K₂₁ in all models regardless of the input SE_y. On the other hand, the entropy values of the correct parameters lg K₁₁ and lg K₁₂ decrease systematically upon lowering the SE_y.

As alternative to the input of a standard deviation to simulate experimental error, Bootstrapping (BS) can be applied. During BS the residuals after fitting are randomly added to the best-fit titration curve to recover the original error. The Shannon entropy for each stability constants after applying Bootstrapping (S = 2000) to each model are visualised in Fig. 5. As observed previously, the Shannon entropy indicates lg K₂₁ as inappropriate stability constant and furthermore the 1:1/1:2 model as the best fitting model.

The application of Monte Carlo simulation after obtaining the best-fit model shows that not only the most probable model should be tested - more information arise from explicit wrong and overfitted models, as in the 2:1/1:1/1:2 sample. The correct model parameter show lower H(x) and σ_dist values, while the incorrect parameters in the complex model are linked to higher values of the descriptors. On the other hand, Monte Carlo Simulation did not indicate that the 1:1 model is not appropriate.

Cross validation

Leave-X-Out Cross Validation (LXO-CV) analysis with X = 1 − 5 was performed and the Shannon entropy and the standard deviation of the histograms obtained for each stability constants were calculated. The Shannon entropy for each LXO-CV run for all stability constants are given in Fig. 6. Individual values obtained per model and parameter follow the trend already obtained in the Monte Carlo simulation. The entropy for the correct parameters (lg K₁₁ and lg K₁₂) tend to be lower than in case of the incorrect stability constant lg K₂₁. Furthermore lg K₁₁ is lowest in the 1:1/1:2 model, followed by the 1:1 model. All this results are obtained with no supplementary input given in contrast to the Monte Carlo Simulation, where a SE_y has to be defined.

The concrete values for H(x) indicate that upon using higher order Cross Validation—leaving more data points out—the entropy and the standard deviation increase. Noteworthy, regardless of the type of the Cross Validation, the statistical descriptors obtained for the model parameters indicate best performance for both lg K₁₁ and lg K₁₂.

Reduction analysis

Reduction analysis was performed on all tested models with the results shown in Fig. 7. Upon continuously removing the very last data point from the data sets in case of the 1:1 model (Fig. 7A), lg K₁₁ constantly increases until a ratio of $\frac{{\left[B\right]}_0}{{\left[A\right]}_0}=1.8$. The best-fit value of lg K₁₁ goes up to 3.2. Afterwards, having only data points below a ratio of 1.8 included, the best-fit value of the stability constants decreases. In contrast to that, the individual parameters in the 1:1/1:2 model do not show such trend above a $\frac{{\left[B\right]}_0}{{\left[A\right]}_0}$ ratio of 2.0. However, using only data points below that ratio, the fitted parameters differ from the originally fitted values. According to that analysis, for correct description of a 1:1/1:2 model at least a ratio of 2.0 is mandatory, which is in line with previous rule-of-thumbs for the Mole-Ratio method. In contrast to the 1:1/1:2 model, the trend in the stability constants of the 2:1/1:1 model show more deviations from the best-fit value upon removing the last data points. The values for lg K₁₁ in the 2:1/1:1/1:2 remain nearly the same above a ratio of 2.0, but differ more below that ratio.

The rationalisation of Reduction Analysis is accomplished using the above introduced partial standard deviation σ_pt (Eq. (11)) where different cut-off values were applied. As cut-off values both a ratio of 1.8 and 2.0 were chosen. The results are represented as bar charts in Fig. 8. Although the values of σ_pt allow rational comparison of reduction analysis, empiricism is mandatory. As the 1:1/1:2 model should recover the simulated data best, the resulting values for σ_pt are again considered as ideal or reference results. Starting with lg K₁₁ and a cut-off ratio of 2.0, the corresponding σ_pt decreases from 0.04037 in the 1:1 model to 0.00724 in the 1:1/1:2 model. In the 2:1/1:1/1:2 model, σ_pt reaches 0.02700 in case of lg K₁₁. In the 2:1/1:1 model, the partial standard deviation for lg K₁₁ exhibits to highest value (0.29115). In analogy, σ_pt for lg K₁₂ increases from 0.03485 about a factor of seven changing the model from 1:1/1:2 to 2:1/1:1/1:2 (0.25491), while σ_pt for lg K₂₁ increases from 0.25578 by a factor of 1.8 to 0.45884. Removing the cut-off, all partial standard deviations get worse and the differences between the models become less clear. Comparing the different sets of σ_pt for both cut-off ratios (1.8 and 2.0) no qualitatively differences can be observed. Slightly lower σ_pt values are obtained for the correct model with the higher cut-off ratio, which is comprehensible since having the lower cut-off ratio, data points are included, at which the saturation concentration for 1:2 complex has not been reached yet. The trend in the partial standard deviations is similar to the trend of Monte Carlo results, where the lowest values of H(x) and σ_dist are obtained for the correct model. Differences occur in the 2:1/1:1 model, where as result of the Monte Carlo simulation both the entropy and σ_dist of the parameter lg K₁₁ are lower than of the parameter lg K₂₁ while using the reduction analysis, the resulting σ_pt is higher for lg K₁₁. Remarkable is the “visual” detection of the insufficient 1:1 model as the removal of data points lead to an increase of the best-fit parameter, each below and above the optimal ratio of 1.0. Even in the 2:1/1:1/1:2 model lg K₁₁ exhibits a more constant behaviour.

Data set with 1:1 stoichiometry

In analogy to the data set based on 1:1/1:2 model, a simulated NMR titration with an underlying 1:1 model was generated and tested using Monte Carlo simulations, Cross Validation and Reduction Analysis. The main results, covering the Shannon entropy of the distribution of individual parameters as well as the partial standard deviation are shown in Fig. 9. A summary of all results is presented in the supplementary material (Figs. S1–S5). The matrix decomposition analysis reveals two factors, which is in agreement with two species contributing to the observed signals (Fig. S20A). Extreme outliers were obtained during both Monte Carlo simulations on top of the best-fit 2:1/1:1 model. Since the Shannon entropy is calculated using 30 bins to partition the histogram and due to the presence of the outliers which enlarge the bins, a reasonable H(x) value for lg K₁₁ was not obtained. The standard deviation is independent of the partition of the histogram, therefore reasonable results can be obtained regardless the distribution. While the combination of the 2:1/1:1 model and Monte Carlo simulations lead to distribution with high standard deviation, the obtained individual parameters after Cross Validation don’t exhibit that pattern. According to the obtained statistical descriptors for each parameter, the 1:1 model describes the simulated titration data best. However, the parameter lg K₁₁, which performs very well in the 1:1 model, performs badly if more complex models are analysed using Monte Carlo simulation or Cross Validation.

Data set with 2:1/1:1 stoichiometry

The collection of the statistical descriptors obtained for analysis of the simulated data set with 2:1/1:1 model are compiled in Fig. 10. A summary of all results is presented in the supplementary material (Figs. S11–S15). The factor analysis indicating two species (Fig. S20C) is not in agreement with the original model using three species contributing to the observed shifts. In contrast to the 1:1 model, where only one parameter was correct, in the current example both, lg K₂₁ and lg K₁₁, are mandatory to describe the data. According to the results visualised in Fig. 10, the descriptors indicating the performance of the parameter of the 2:1/1:1 model are better (lower) than in the 2:1/1:1/1:2 model. Furthermore the descriptors indicate that lg K₁₁ performs in the 2:1/1:1 model better than in the 1:1/1:2 model. However, in contrast to the previously discussed data of a simulated 1:1/1:2 experiment, the values of the descriptors are in general higher, although the absolute values are of limited meaning. The differences of the descriptors for lg K₁₁ in either the applied 1:1 or the 2:1/1:1 model are not as clear as between the applied 1:1 and the applied 1:1/1:2 model. In case of H(x) for lg K₁₁ after MC with 2000 steps, the entropy drops from 0.08374 (1:1 model) to 0.05235 (1:1/1:2 model) by a factor 0.63 in the simulated 1:1/1:2 titration and from 0.25172 (1:1 model) to 0.22620 (2:1/1:1 model) by a factor of 0.90 in the simulated 2:1/1:1 titration data set. Using the L3O-CV results, the difference is clearer. The entropy drops from 0.50499 (1:1 model) to 0.17749 (2:1/1:1 model) by factor 0.35 in case of the 2:1/1:1 titration. The change of the entropy in the 1:1/1:2 titration is by a factor of 0.68, from 0.04696 (1:1 model) to 0.03214 (1:1/1:2 model) respectively. However, the results in Fig. 10 indicate that the 1:1/1:2 model performs worse than the 2:1/1:1 or 1:1 model, therefore only 1:1 model could be considered as alternative model. Comparing this results with results obtained for the pure 1:1 model, lg K₁₁ performs clearly worse if a 2:1/1:1 model is fitted to a 1:1 data set. Following this pattern, the statistical post-process did not fail to detect 2:1/1:1 stoichiometry correctly, although the results are not as clear as in case of the 1:1/1:2 model.

Data set with 2:1/1:1/1:2 stoichiometry

In contrast to the correct detection of the appropriate models in the previous discussion, the individual parameters obtained for the stability constants don’t reflect lg K₂₁ as suitable parameter (Fig. 11). A summary of all results is presented in the supplementary material (Figs. S16–S20). According to the factor analysis, three species contribute to the observed signals (Fig. S20D), which is not in agreement with the data set. In both models, 2:1/1:1 and 2:1/1:1/1:2, the statistical descriptors for lg K₂₁ show higher values than the descriptors of lg K₁₁ and lg K₁₂. Although the identification of the correct stoichiometric model was not possible with the statistical post-processing, this can be rooted back to the simulated data. The titration data was simulated with lgK₂₁ = 1.85, lgK₁₁ = 3.60 and lgK₁₂ = 2.40, where only the stability constants lg K₁₁ and lg K₁₂ were recovered correctly (3.61 and 2.42). The stability constants lg K₂₁ was estimated to be 0.21, indicating no 2:1 stoichiometry at all. The correct estimation of 2:1 stoichiometry is somewhat difficult, as the concentration of the 2:1 species is lower than for the other species. Therefore the influence of this complex to the overall observed signal is of less importance.

Averaged statistical descriptors

Averaged statistical descriptors for the stability constants calculated for each model and for each parameter are illustrated in Figs. 12 and 13. The chart in Fig. 12 shows the average of the descriptors like H(x) and σ_dist obtained for all parameters within a single model. For example the blue bar represents the average of the descriptors for all parameters of the 1:1 model, which is lg K₁₁. The orange bars represent the average of the descriptors for lg K₁₁ and lg K₁₂, which are the stability constants in case of the 1:1/1:2 model. As each bar represents the average values of the descriptors for each applied model, the general quality of the fit using a special model on top of the data set is easily deduced. The model fitted best to the simulated 1:1 systems is the original 1:1 model, where however the average descriptor is obtained for only one stability constant. The average partial standard deviation in case of the other tested models is higher, with the highest value obtained for the 1:1/1:2 model. According to the average σ_pt, the 1:1/1:2 system is best described by the 1:1/1:2 model and worst by the 2:1/1:1 model. A similar agreement between the true simulated model and lowest averaged descriptors was obtained in case of the 2:1/1:1 model, however the difference of the σ_pt value compared to the 1:1 model is very small. Using the average descriptors for σ_pt, the 2:1/1:1/1:2 model was not identified as the correct one. Better results were obtained fitting a 1:1 and 1:1/1:2 model to the simulated data with 2:1/1:1/1:2 stoichiometry, which is in accordance with results discussed in Section ‘Data set with 2:1/1:1/1:2 stoichiometry’. Results obtained with Reduction Analysis were comparably recovered using Cross Validation and calculated Shannon entropy and standard deviation of the distribution. Comparing the averaged σ_pt values with the averaged Shannon entropy, the difference is similarly small in case of the 1:1 and 2:1/1:1 model applied to the 2:1/1:1 data set (Figs. 12A and 12B). In contrast, using σ_dist to characterise the distributions of stability constants the standard deviation obtained for both, the 1:1 and 2:1/1:1 model differ more clearly (Fig. 12C). The results with the remaining CV and MC calculations can be found in the supporting information. Apart from the mixed 2:1/1:1/1:2 model, post-processing after RA, CV and MC identifies the correct models by assigning them the lowest value of the averaged descriptors compared to the remaining models. The absolute values of these descriptors have not been identified as meaningful on its own.

Similar to the average descriptors obtained for each model individually, the average can be calculated for the all parameters, e.g., stability constants. In Fig. 13A, the blue bar represents the average of four σ_pt values obtained for the stability constant lg K₁₁ from the 1:1, the 2:1/1:1, the 1:1/1:2 and the 2:1/1:1/1:2 model. The green bar represents the average of all lg K₂₁ stability constants, however obtained for fewer models. The remaining average descriptors are obtained in a similar fashion for all stability constants of all data sets. It is obvious from the charts in Fig. 13 that apart from the 2:1/1:1/1:2 model, the most appropriate stability constants in general have lower average values compared to those which are inappropriate. This is true for RA as well as L3O-CV using both descriptors. While this recovers the expectation derived by the previous discussion, there is a fundamental shortcoming in this context. In contrast to the average descriptors for each model, focusing on the stability constants or the other model parameters, an a priori knowledge of the number of model parameters is mandatory. Given the result in Fig. 13A for the 1:1 model, the difference between σ_pt for lg K₁₁ and lg K₂₁ is much lower than compared to σ_pt for lg K₁₂. On the other hand, in case of the 1:1/1:2 model, the difference between lg K₁₂ and lg K₂₁ is much lower than in the simulated 1:1 model. And the difference between the σ_pt values for the correct stability constants lg K₁₁ and lg K₁₂ (1:1/1:2 model) is comparable to the difference between the correct lg K₁₁ and the inappropriate lg K₂₁(1:1 model). Some differences become clearer if the CV results were analysed. However, if both, the number of appropriate parameters and the parameter itself, are unknown, the parameter-wise average of the descriptors is not meaningful.

Based on the results for both strategies to calculated average descriptors, the meaningful approach is to average over all parameters for a given single model and not to average over all identical parameters for all models. However, both strategies are applied automatically upon comparing the results of the statistical post-processing for applied models in SupraFit.