Separating the effects of experimental noise from inherent system variability in voltammetry: the [Fe(CN) 6 ] 3−/4− process
Recently, we have introduced the use of techniques drawn from Bayesian statistics to recover kinetic and thermodynamic parameters from voltammetric data, and were able to show that the technique of large amplitude ac voltammetry yielded significantly more accurate parameter values than the equivalent dc approach. In this paper we build on this work to show that this approach allows us, for the first time, to separate the effects of random experimental noise and inherent system variability in voltammetric experiments. We analyse ten repeated experimental data sets for the [Fe(CN) 6 ] 3−/4− process, again using large-amplitude ac cyclic voltammetry. In each of the ten cases we are able to obtain an extremely good fit to the experimental data and obtain very narrow distributions of the recovered parameters governing both the faradaic (the reversible formal faradaic potential, E 0 , the standard heterogeneous charge transfer rate 1 constant k 0 , and the charge transfer coefficient α) and non-faradaic terms (uncompensated resistance, R u , and double layer capacitance, C dl ). We then employ hierarchical Bayesian methods to recover the underlying "hyperdistribution" of the faradaic and non-faradaic parameters, showing that in general the variation between the experimental data sets is significantly greater than suggested by individual experiments, except for α where the inter-experiment variation was relatively minor. Correlations between pairs of parameters are provided, and for example, reveal a weak link between k 0 and C dl (surface activity of a glassy carbon electrode surface). Finally, we discuss the implications of our findings for voltammetric experiments more generally.
Introduction<!>Experimental methods<!>Mathematical Modelling<!>Parameter Recovery Methods<!>Bayesian Inference<!>Calculating the likelihood<!>Markov Chain Monte Carlo parameter inference<!>Practical implementation of the Metropolis-Hastings Algorithm<!>Hierarchical Bayesian Inference<!>Generating synthetic data as a test case<!>Synthetic data<!>Experimental data<!>Discussion<!>Conclusion
<p>In a previous paper, 1 we described the use of Bayesian inference for quantitative comparison of voltammetric methods for investigating electrode kinetics. We illustrated the utility of the approach by comparing the information content in both dc and ac voltammetry at a planar electrode for the case of a quasi-reversible one electron reaction mechanism. Using both synthetic and experimental data, we were able to demonstrate that realistic levels of purely random experimental (Gaussian) noise have a relatively minor affect on the inverse problem of recovering both the faradaic (the reversible formal potential, E 0 , the standard heterogeneous charge transfer rate constant k 0 , and the charge transfer coefficient α) and non-faradaic (uncompensated resistance, R u , and double layer capacitance, C dl ) parameters that govern this reaction mechanism. We also demonstrated the clear advantages in terms of accuracy of parameter recovery of the large amplitude ac approach. With this ability of being able to recover parameter values reliably from a single experimental data set in place, we are now in a position to go on to investigate system level variability i.e. if we repeat the same experiment multiple times and implement our Bayesian parameter recovery procedure for each data set independently, how reproducibly do we recover the governing parameters?</p><p>To approach this problem, we again return to the "pathological" [Fe(CN) 6 ] 3−/4− process.</p><p>It is well known that for this system the kinetic parameters reported are highly variable even when using apparently identical electrodes and conditions (see [2][3][4][5][6][7][8][9][10][11][12][13][14] and references cited therein). Using our new approach we are able to show that this difficulty is not due to the impact of random experimental noise, since for each individual data set we are able to fit the mathematical model to the experimental data extremely accurately. However, we are able to demonstrate that the recovered values from each individual data set vary significantly from one another, i.e. the system itself varies between experimental runs, and so the recovered parameter values are extremely sensitive to the precise experimental conditions pertaining on that particular run. Since we have already demonstrated its advantages over the dc approach, we restrict ourselves to the ac case in this paper.</p><!><p>Details of the experimental data sets used in this paper have been given previously in. 2,15 In summary, large amplitude ac voltammetry was performed in a standard three-electrode cell, using a glassy carbon macrodisk (diameter 3 mm) working electrode. All potentials are reported versus an Ag/AgCl/KCl(3 M) reference electrode (hereinafter Ag/AgCl). The frequencies, amplitude and scan rate were 9.02 Hz, 0.080 V and 0.894 Vs −1 respectively and data were collected over the potential range of 0.5 to -0.1 V versus Ag/AgCl. The surface area of the electrode was estimated as 0.070 cm 2 . The value of the diffusion coefficient, D, of [Fe(CN) 6 ] 3− was found to be 7.2 × 10 −6 cm 2 s −1 , as described in. 15 In this paper we make use of the ten repeated data sets for ac voltammetry taken from 15 for the reduction of aqueous 1.0 mM [Fe(CN) 6 ] 3− in 3 M KCl aqueous electrolyte solution. In order to reduce the amount of computation required for the parameter inference, we use a moving average window of length 21 to reduce the number of experimental data points within each data set to about 25,000 data points.</p><!><p>The details of the mathematical and computational modelling approach that we have taken in this paper were given previously in. 1,15 In summary, our chosen experimental system is modelled as a quasi-reversible reaction</p><p>where species A and B are in solution, and E 0 , k 0 , and α are the reversible formal potential, standard heterogeneous charge transfer rate constant at E 0 and the charge transfer coefficient, respectively. We assume that the Butler-Volmer formalism applies to the electron transfer process. [16][17][18] We also assume that both convection and migration can be neglected since we are using a macrodisk stationary electrode and an excess of supporting electrolyte, respectively. Since we assume equal diffusion coefficients for each species A and</p><p>where x is distance from the electrode surface and t is time. The initial and boundary conditions are</p><p>At the electrode surface, x = 0, for t > 0, we have the conservation and flux conditions</p><p>along with the Butler-Volmer condition</p><p>Here, I f is the faradaic current, S is the electrode area, and E eff (t) is the effective applied potential (defined below).</p><p>We complete the model by defining E app (t) to be the applied potential, then for the case of an ac voltammetry ramp we have</p><p>where v is the sweep rate, E start is the initial potential, t reverse is the time of switching from the forward to the reverse sweep in cyclic voltammetry, ω is the radial frequency and ∆E is the amplitude of the sine wave. The effective applied potential can now be defined as</p><p>where E drop models the effect of uncompensated resistance, R u . I tot is the total (measured) current, and combines the faradaic current and the background capacitive current, I c , which can be modelled as</p><p>where C dl is the double layer capacitance (assumed constant in this work), and then</p><p>Equations 2-8 are non-dimensionalised as described previously. 1,15 The resulting nondimensional system of equations is solved using an implicit finite difference method with an exponentially expanding grid (again, as described previously 19 ). We can now see mathematically that the reaction mechanism in Eq. 1 is governed by five parameters (E 0 , k 0 , α, C dl , R u ), and we will collectively denote these parameters by the vector θ. The inverse problem that we wish to solve can be defined as finding the best possible approximation to θ given our measured experimental output trace of the current I data tot versus potential.</p><!><p>In a recent paper we described in detail how methods of Bayesian inference can be used to solve the inverse problem of parameter recovery 1 from voltammetric data. We illustrated how these methods yield not only a point estimate for each parameter of interest, but also a measure of our confidence in that estimate. Full details of the approach that we adopted, including the algorithms used can be found in, 1 so that here we simply give a brief outline of the methods.</p><!><p>We denote by y = (y 1 , . . . , y T ) the experimental data trace, that is, the total measured current I data tot, t at each time point t. In our previous work (see Figure 4 and Table S2 of 15 ) we demonstrated that, to a very good approximation, experimental measurements can be assumed to be subject to normally distributed random noise, typically with zero mean and some standard deviation which we will denote by σ (we showed that typical values of the standard deviation of the experimental noise are in the range 1 to 2% of the peak current). Our mathematical model of the system, Eqs 2 to 8 above, assumes that the observed experimental data y is a function of the parameters of interest θ, θ = (E 0 , k 0 , α, C dl , R u ).</p><p>We now assume further that the parameters θ are themselves drawn from a probability distribution. We can then frame our inverse problem as trying to find this probability distribution for θ given the observed values of the data y, and denote this probability distribution as P (θ|y) (the vertical line indicates that the values of y are given). In Bayesian inference, P (θ|y) is termed the posterior probability density or posterior distribution; this is the distribution that we want to approximate.</p><p>We now make use of Bayes' rule which states P (θ|y) = P (y|θ)P (θ)</p><p>where P (θ) is called the prior distribution of θ and is chosen to capture any prior knowledge we have of θ before any experimental observation. The distribution P (y|θ) is the probability density of the experimental data y given a model parameterised with parameters θ, and is termed the likelihood of the data; assuming a known distribution of the error in the data this likelihood can be calculated. P (y) is a normalising term (which is the integral of all possible densities P (y, θ) = P (y|θ)P (θ) over all values of θ), and ensures that the posterior density P (θ|y) integrates to 1. In practice, the calculation of this normalising term (which can be very computationally expensive) is avoided by considering ratios of the likelihood (see below).</p><!><p>Writing the likelihood as</p><p>we can re-arrange Bayes' rule in Eq. 9 to give</p><p>Since we assume that the errors are independent at each time point, the conditional probability density of observing the experimental trace from t = 1, . . . , T given θ is simply the product of the probability density functions at each time point, that is, the likelihood is given by</p><p>using the assumption that the experimental noise is normally distributed with a mean zero and variance of σ 2 . For notational simplicity we have set</p><p>Since, in most experimental situations in electrochemistry, a priori we will have only a rough idea of what the values of the parameters are likely to be, we assume an "uninformative" prior and use a uniform distribution for each parameter across a suitably wide range, that is</p><p>where c is a non-zero finite normalizing constant. The bounds for this hypercube were set to</p><p>Note that this prior is only used for the single level Markov Chain Monte Carlo (MCMC) algorithm; in the hierarchical MCMC algorithm, this is replaced by a multivariate normal, as described below. However, these bounds are still used in the hierarchical MCMC algorithm to prevent the lower level samplers from accepting any samples that lie outside these bounds.</p><!><p>To obtain a sample from the posterior distribution P (θ|y) we make use of the Markov Chain Monte Carlo method. In outline, this involves finding an approximation to the posterior distribution P (θ|y) by drawing a finite (but sufficiently large to be accurate) number of samples from this distribution. To do this we simulate a Markov Chain whose limiting distribution is the required posterior distribution using an efficient implementation of the Metropolis-Hastings algorithm, 20 within which candidate parameter sets are proposed from a proposal distribution q(θ cand |θ i ) which depends only on the previously accepted parameter set θ i ;</p><p>we take q(θ cand |θ i ) to be a multivariate normal distribution. If the proposed parameter set contains any parameters outside the range of the prior, then the parameter set is assigned an acceptance probability of 0, i.e. it is rejected, and the previously accepted parameter set is added to the Markov chain -that is, θ i+1 = θ i . Otherwise, we compare θ cand to the current parameter set θ i by calculating the ratio of the posteriors of the two parameter sets. If the candidate parameter set has a greater posterior density value than the existing parameter set then it will be added to the Markov chain, that is θ i+1 = θ cand . Otherwise, (making use of Eq. 11) the proposed parameter is accepted with probability, r, given by</p><p>If the proposed parameter set is rejected (with probability 1 − r), then the previously accepted parameter set is again added to the Markov chain -that is, θ i+1 = θ i .</p><p>Since the likelihoods for large samples are extremely small, in practice we work with the natural log of the likelihood which reduces to</p><p>where terms which are constant in θ have been removed (since these will cancel on taking the difference of log-likelihoods in the Metropolis-Hastings algorithm). More comprehensive descriptions of the theory of MCMC can be found in the statistics literature (see for example 21 ).</p><!><p>In practice, we make use of an adaptive covariance matrix version of the Metropolis-Hastings algorithm which helps identify the directions in parameter space which have the highest likelihood values. 22 At each iteration of the algorithm, the covariance matrix of the multivariate normal distribution is updated and a scalar value is also updated to define the width of the distribution. In the results presented in this paper, we run our MCMC chains for 10,000 samples and discard the first 5,000 samples as 'burn in' (see 21 ). To ensure efficiency in our Monte-Carlo sampling, we first find the location of the optimum, that is, the maximum likelihood estimate of θ using a standard global minimisation algorithm (we use the cma-es algorithm 23 ), which we use as a seed point for the MCMC algorithm as described in. 1 In the results section, these samples are shown as histograms which illustrate the nature of the posterior distribution.</p><!><p>In the results section, the MCMC algorithm described above is used to estimate the parameter values and their posterior distributions from each of the ten repeat runs of the ac voltammetry experiment for the reduction of aqueous 1 mM [Fe(CN) 6 ] 3− , as described in the Experimental methods section. This will allow us to show the degree of variability in the recovered values of the parameters across these ten data sets. This will in turn allow us to postulate that on each repeat of the experiment the parameters θ = {E 0 , k 0 , α, C dl , R u } are themselves drawn from a multivariate normal distribution with mean hyper-parameters µ = ( Ê0 , k0 , α, Ĉdl , Ru ) and a covariance hyper-parameter matrix Σ. Our aim now is to sample from the distributions of µ and Σ to enable us to quantify how θ varies between different experiments.</p><p>We first write the posterior distribution of all of our parameters given the data, taking into account our new hierarchical model structure</p><p>= P (µ, Σ)</p><p>where θ 1 , ..., θ n are all the bottom level parameters for the n different experiments, and y 1 , ..., y n are the corresponding measurements.</p><p>We choose a multivariate normal distribution for the bottom level parameters θ</p><p>It now remains to choose a suitable hyper-prior P (µ, Σ), which for ease of computation is generally taken to be a normal-inverse-Wishart distribution (see for example 24 )</p><p>The distributions in Eqs. 19 and 20 are chosen to be conjugate, so that the conditional distribution of the hyper-parameters can be analytically derived as 25</p><p>where θ and C are the sample mean and covariance of the bottom level parameters</p><p>parameters, combined with classical Gibbs sampling and Eq. 21 to sample from the hyperparameters (see Algorithm 1 for details). For the hyperprior parameters we use κ 0 = 0 and ν 0 = 1. We set µ 0 to be the centre point of the 5-dimensional hypercube in Eq. 14, and Ψ as a diagonal matrix with the standard deviation of each parameter set to half the width of this same hypercube. To prevent the bottom level MCMC chains from exploring unphysical parameter regimes, we automatically reject any proposed point which lies outside of the hypercube.</p><p>Algorithm 1 Metropolis within Gibbs s = 0 µ s , Σ s = SampleFromN IW(µ 0 , κ 0 , ν 0 , Ψ) {using Eq. 20} repeat for i = 1 to n do θ i = AdaptiveMCMCStep(µ s , Σ s ) {see Eq. 15 and enclosing section} end for</p><!><p>To test our inference procedure and algorithms we make use of "synthetic" test data i.e. we generate ten synthetic data sets by solving Equations 2 to 8 for a randomly chosen set of values of θ = (E 0 , k 0 , α, C dl , R u ), drawn from a multivariate normal distribution with mean θ true = (7.27, 2.01, 0.53, 3.70 × 10 −3 , 1.06 × 10 −2 ) (in non-dimensional units), and standard deviation σ true = (0.06, 0.7, 0.005, 0.7 × 10 −3 , 0.3 × 10 −2 ). We then add randomly generated Gaussian noise at each time point with zero mean and standard deviation of 0.3% of the maximum current (chosen to match the experimental noise level).</p><p>13</p><!><p>To verify that we have implemented our hierarchical Bayes algorithms correctly we briefly describe the results of using synthetic data to test our inference procedure. Figure 1 shows the sampled distributions of the hyper-parameters. The dashed lines indicate the expected peak of each distribution, as calculated from the sample mean and variance of the ten true values of θ true . As can be seen the hierarchical MCMC algorithm obtains the correct peak for all the five different mean and variance hyper-parameters.</p><!><p>Figure 2 shows an example of the initial fitting process for experimental data set 1. Once the best fit for each data set was obtained using maximum likelihood estimation and the cma-es optimisation algorithm, each of the bottom level MCMC algorithms was initialised at these points. Then Algorithm 1 was used to generate 5,000 samples (10,000 total samples, of which the first 5,000 was discarded as burn-in) of the parameters θ i and hyper-parameters µ and Σ.  ). Simulation parameters were obtained by maximum likelihood estimate, and their values in dimensional units were E 0 = 0.214 V, k 0 = 0.010 cm s −1 , α = 0.528, C dl = 16.9 µF cm −2 , R u = 0.00 Ω.</p><p>taken using the original non-hierarchical model (i.e. just running ten independent MCMC chains on the ten different data sets), and this is shown in Figure 3 (right column). As can be seen, the result in this case are almost identical to the hierarchical model, giving us confidence that we are capturing the distributions of θ i correctly in each case.</p><p>The chief benefit of the hierarchical model is that it allows us to quantify (with the hyper-parameters) the variability of the parameters between different experiments. Figure 4 shows the histograms for the hyper-parameters samples of µ (left) and Σ (right). Also shown as vertical dashed lines are the the sample mean and variance of the concatenated ten θ i chains. The most obvious feature of these plots is that these hyper-parameter distributions clearly show a much wider confidence interval as compared with the individual θ i distributions. While the value of each parameter for each individual experiment is known with high accuracy, once the variability between experiments is taken into account we see that the possible parameter range is significantly broadened.</p><!><p>The Bayesian data analysis approach introduced in this paper provides access to fundamentally new knowledge that assists in elucidating nuances that have contributed to the highly non-reproducible electrode kinetic data published for the "pathologically variable"  The left column plots show the samples obtained using the hierarchical model, and these also show the posterior predictive distribution P (•) for each variable, calculated by summing the Gaussian distributions described by the samples of the hyper-parameters µ and Σ (see Figure 4 for histograms of these hyper-parameter samples). For comparison, the right column plots show the samples of θ i with no hierarchical model (and thus no hyper-parameter samples). This is identical to the analysis done in our previous paper.   0.002 to 0.018 cm s −1 is evident in 10 individual experiments at a nominally identical electrode surface (Figure 3), conformance to the quasi-reversible model is exceptionally good for each individual experiment.</p><p>The variation in performance of the now very widely used glassy carbon electrode was identified as a point of concern soon after the material was introduced into electroanalytical chemistry. 26 In essence, the exact nature of the glassy carbon, and indeed other carbon based surfaces, has their origin in the method (e.g. temperature) of manufacture, nature of pre-treatment and its history as an electrode (see for example 13,14,[26][27][28][29][30][31][32] and references cited therein). Accordingly, it seems likely that the widely variable electrode kinetic data associated with published studies on the [Fe(CN) 6 ] 3−/4− process (summarised in Eq. 24)</p><p>mimics the variability of the surface state used in the different publications.</p><p>In practice, carbon electrodes are highly heterogeneous with surface defects and organic functional groups in abundance when used in aqueous electrolyte media. [28][29][30] This seems to translate into an electron transfer process occurring at a surface consisting of microscopically small and distinctly different regions that must be sufficiently close so that complete overlap of diffusion layers occurs on the measurement timescale, 33 presumably allowing the entire surface to be successfully modelled by approximating the mass transfer by planar diffusion.</p><p>Thus, even though highly heterogeneous with variable k 0 values at the microscopic level, the surface behaves as though it is fully homogeneous, hence giving rise to measurement of an apparently single "averaged" k 0 value.</p><p>Electron transfer reaction mechanisms are often classified into inner and outer sphere categories in both homogeneous chemical redox reactions that occur in the solution phase and in heterogeneous reactions in electrochemistry that take place across an electrode-solution interface. 16,33 In a chemical redox reaction, the homogeneous outer sphere class mechanism * Fe(CN)</p><p>which is the homogeneous analogue of the electrochemical one of interest in this study, while in the inner class, ligands involving bridging to a common metal centre may be involved.</p><p>In an outer sphere electrochemical process, the plane of closest approach to the Outer Helmholtz Plane does not allow penetration of the layer of non-specifically adsorbed or coordinated solvent adhered to the electrode surface by reactants. In contrast, inner sphere electron transfer processes involve specifically adsorbed reactants, 16,33 and therefore are anticipated to exhibit electrode kinetics that are strongly dependent on the chemical nature of Before attempting to identify the factors that may be most relevant to the variation encountered in the electrode kinetics, it is worthwhile reviewing the electrode pre-treatment regime. A nominally 3 mm diameter GC working electrode was purchased from BAS. Prior to each set of measurements the surface of the GC electrode was thoroughly polished with 0.3 µm alumina powder in the form of an aqueous slurry on a wet polishing cloth (BAS).</p><p>After polishing, the electrode was repeatedly washed with high purity water and subjected to sonication for 10 to 20 s. To provide effective removal of any residual alumina powder, after the initial sonication, the electrode was carefully wiped with a clean wet polishing cloth, again washed with water and sonicated for 10 to 20 s in high purity water. Finally, the working electrode was dried under a nitrogen stream. As far as possible, each electrode preparation was undertaken under identical conditions; nevertheless, the electrode kinetics differ substantially from experiment to experiment. However, of course the history of the electrode could be important since the electrode used for experiment ten will have a more extensively polished electrode than that used in experiment one.</p><p>The narrow distribution of each P i (•) in Figure 3 enables us to conclude that the contribution to parameter variability from noise associated with each particular experiment is minimal for all five parameters estimated. From perusal of Figures 3 and 4, we can conclude that while the degree of experiment-to-experiment variation changes with respect to each parameter, it is in general significantly greater than the parameter variability from noise alone (i.e. the width of each P i (•)). From the hyper-parameter samples we can calculate the distribution P (•) of each parameter that would be expected from another repeat experiment (shown in Figure 3, left column). Quantitatively, the parameter mean values (with one standard deviation values in parenthesis) as deduced from this analysis are as follows:</p><p>E 0 = 0.215 V vs Ag/AgCl (0.002), k 0 = 0.010 cm s −1 (0.005) α = 0.521 (0.006), C dl = 20.1 µF cm −2 (4.5), R u = 15.2 Ω (21.0). Intriguingly, experiment-to-experiment variation in α is very small. That is, its value lies in a very narrow range of about 0.515 to 0.525 with a standard deviation of just 0.006. It is also notable that while the value of R u covers a range from zero to about 35 Ω, it is always small. This means that the more important ohmic drop term that can distort voltammograms also is small and hence not highly significant.</p><p>The large standard deviation of 21.0 Ω with a mean value of 15.2 Ω allows us to conclude the R u parameter only has a minor impact on the ac voltammetry (highly conducting 3 M KCl aqueous electrolyte) and hence on E 0 , k 0 , or α variability. Although the experimentto-experiment variability of E 0 is much greater than predicted by P i (E 0 ), it is still fitted to a tight regime (about 2 mV) when compared with the total voltage scan range. Indeed, in the absence of adsorption, E 0 is theoretically predicted to be completely independent of electrode material or state, suggesting that the inter-experiment variability of E 0 is due to compensating factors from other parameters. C dl displays some variability with values lying within the range of about 16 to 23 µF cm −2 which confirms that the electrode is not in an identical surface state, and hence does have a variable level of activity for each experiment.</p><p>k 0 as noted above encompasses a wide range of about 0.002 to 0.018 cm s −1 (Figure 3) with a mean value of 0.010 cm s −1 and a standard deviation of 0.005 cm s −1 . The question that arises is whether there is a correlation between C dl or surface activity and k 0 .</p><p>Figure 5 shows quantitatively the correlation between pairs of parameters using scatterplots of the mean hyper-parameter µ samples. Clearly, R u and α are not correlated at all with each other or any other parameter (i.e. they give "shot gun" correlation plots).</p><p>However, there appears to be a weak correlation of k 0 and C dl with larger k 0 values being associated with data sets having higher C dl values. There also appears to be a weak correlation between k 0 and E 0 , with the more positive E 0 values coinciding with the larger C dl values. One needs to be careful not to over interpret the significance of weak correlations, as other non-quantified parameters also can be operative. However, a larger capacitance current caused by a larger activation of the glassy carbon surface (more surface functionality) may be expected to increase the rate of an inner sphere process, like the [Fe(CN) 6 ] 3−/4− one probed in this work. The origin of the weak correlation of k 0 and E 0 is more problematical.</p><p>The presence of extremely weak adsorption not accommodated in the model is one possible explanation, but a small drift in reference electrode potential cannot be ruled out.</p><!><p>In summary, the Bayesian inference-inspired strategy introduced in this paper for data evaluation represents a significant advance in understanding the contribution of different parameters to voltammetric data and their significance in experiment-to-experiment variability at a heterogeneous electrode in a manner that has not been possible in earlier studies. Perhaps transfer and mass transport by planar diffusion even though the variation in k 0 is quite substantial. Intriguingly, α, unlike k 0 , does not reveal significant experiment-to-experiment variation in this data analysis exercise. This may be expected if conformance to the Butler-Volmer model is strong.</p><p>In early electrode kinetic studies, the ideal and very homogenous mercury electrode was used. It is now evident that electrode design is becoming very sophisticated, particularly when advances in materials science are used to generate highly heterogeneous electrode materials. 34 Thus, there is a tendency in electrochemistry nowadays to use far more complex electrodes than the historically important liquid mercury and pure solid metal surfaces. The new breed of highly heterogeneous electrodes imply that data analysis strategies in the future will also need to be more sophisticated. That is, use of models significantly more complex than used in this study with Butler-Volmer theory for electron transfer and mass transport by planar diffusion, as well as Bayesian forms of data analysis, will become increasingly essential if informative experimental versus simulation comparisons are to be reported. Distributions for all ten lower level adaptive MCMC samples of the reversible formal potential E 0 . The coloured histograms show the lower level posterior for each of the 10 repeat experiments P i (E 0 ). The solid line shows the predictive posterior for E 0 , which describes the distribution of each parameter that would be expected from another repeat of the experiment.</p>