Article pubs.acs.org/ac
A Comparison of Fully Automated Methods of Data Analysis and Computer Assisted Heuristic Methods in an Electrode Kinetic Study of the Pathologically Variable [Fe(CN)6]3−/4− Process by AC Voltammetry Graham P. Morris,† Alexandr N. Simonov,‡ Elena A. Mashkina,‡ Rafel Bordas,§ Kathryn Gillow,† Ruth E. Baker,† David J. Gavaghan,*,§ and Alan M. Bond*,‡ †
Mathematical Institute, University of Oxford, Andrew Wiles Building, Radcliffe Observatory Quarter, Woodstock Road, Oxford, OX2 6GG, United Kingdom ‡ School of Chemistry, Monash University, Clayton, Vic. 3800, Australia § Department of Computer Science, University of Oxford, Wolfson Building, Parks Road, Oxford, OX1 3QD, United Kingdom S Supporting Information *
ABSTRACT: Fully automated and computer assisted heuristic data analysis approaches have been applied to a series of AC voltammetric experiments undertaken on the [Fe(CN)6]3−/4− process at a glassy carbon electrode in 3 M KCl aqueous electrolyte. The recovered parameters in all forms of data analysis encompass E0 (reversible potential), k0 (heterogeneous charge transfer rate constant at E0), α (charge transfer coefficient), Ru (uncompensated resistance), and Cdl (double layer capacitance). The automated method of analysis employed time domain optimization and Bayesian statistics. This and all other methods assumed the Butler−Volmer model applies for electron transfer kinetics, planar diffusion for mass transport, Ohm’s Law for Ru, and a potential-independent Cdl model. Heuristic approaches utilize combinations of Fourier Transform filtering, sensitivity analysis, and simplex-based forms of optimization applied to resolved AC harmonics and rely on experimenter experience to assist in experiment−theory comparisons. Remarkable consistency of parameter evaluation was achieved, although the fully automated time domain method provided consistently higher α values than those based on frequency domain data analysis. The origin of this difference is that the implemented fully automated method requires a perfect model for the double layer capacitance. In contrast, the importance of imperfections in the double layer model is minimized when analysis is performed in the frequency domain. Substantial variation in k0 values was found by analysis of the 10 data sets for this highly surface-sensitive pathologically variable [Fe(CN)6]3−/4− process, but remarkably, all fit the quasi-reversible model satisfactorily.
T
DC cyclic voltammetry is probably the most widely used method to study electrode kinetics. It is excellent for initial qualitative assessment but not straightforward in quantitative studies. Many experiments over a wide range of scan rates are needed, and maintaining a reproducible electrode surface for the duration of the experiment is problematic. A very sensitive technique for the quantitative evaluation of electrode kinetics is based on the AC version of voltammetry in which a periodic signal is superimposed onto a DC potential ramp of the kind used in DC cyclic voltammetry. In recent work from these laboratories,6−9 when using large amplitudes that give rise to substantial nonlinearity in the response, it has been shown that
he technique of voltammetry is used widely in many branches of chemistry and other disciplines in order to study electron transfer processes, with the reactant and product being most commonly soluble in the solution phase. Ideally, theoretical and experimental comparisons would be undertaken routinely so that a quantitative account of the significance of the measurements undertaken is provided. For solution phase studies, the theory for voltammetry is now very well developed, and user-friendly simulation packages are available that encompass a large range of mechanisms, e.g., DigiSim,1 DigiElch,2 COMSOL Multiphysics,3 KISSA,4 and MECSim.5 Nevertheless, quantitative reports of electrode kinetics and the related parameters are still not routine, and even when undertaken it is still unusual to find a detailed analysis of the quality of agreement between theory and experiment. © 2013 American Chemical Society
Received: July 18, 2013 Accepted: October 27, 2013 Published: October 28, 2013 11780
dx.doi.org/10.1021/ac4022105 | Anal. Chem. 2013, 85, 11780−11787
Analytical Chemistry
Article
[Fe(CN)6 ]3 − + e− ⇌ [Fe(CN)6 ]4 −
analysis of the current−time data can be aided by use of a Fourier transform−inverse Fourier transform sequence of data analysis. Fourier transformation converts experimentally obtained time domain data to the frequency domain. When data are displayed as a power spectrum, band selection can be undertaken on the resolved DC and AC components. Inverse Fourier transformation can then be used to convert the resolved components back into the time domain and the comparisons of simulated and experimental responses undertaken on each of the aperiodic DC and AC harmonic components. Only a single experiment is needed in this approach, so the electrode surface is invariant with respect to impact on data analysis. The simplest electrode process available for analysis can be written as A ⇌ B + e−
(2)
which is very electrode-surface dependent,10−17 are considered. On carbon surfaces, the reaction in eq 2 can be regarded as pathologically variable in the sense that reproducibility of literature outcomes apparently using the same electrode pretreatment has been almost impossible to achieve in different laboratories.11−13,15,17 An extensive reappraisal of the voltammetry of the [Fe(CN)6]3−/4− process at carbon electrodes has been published recently by Unwin et al.,17 who have scrutinized the relationship between the history of the electrode surface and the electrode kinetics at the macro-, micro-, and nanoscale levels. Other papers have also examined the problems of reproducibility at carbon surfaces.18−20 For further details on the origin of the pronounced dependence of the kinetics of reaction 2 on the characteristics of the electrode, the reader is referred to these papers and references cited therein. In essence, the surface of an electrode constructed from a material like glassy carbon is kinetically heterogeneous, with surface functional groups and defects being crucial to the end result. Abundant literature exists on how to interpret data for the domains of behavior at carbon electrodes.21,22 However, as it emerges from the FTACV studies in this paper, the theoretical treatment of the [Fe(CN)6]3−/4− process as quasi-reversible is appropriate, but reproducibility of k0 values detected from the analysis of a series of experiments conducted under apparently identical conditions is ultimately poor. For the first time, in this paper direct analysis of raw data is considered as an alternative to analysis of the individual DC and AC harmonics resolved by Fourier transform methods. Additionally, the role of noise and reproducibility of data are considered, and conclusions related to preferred methods of data presentation and analysis are drawn. Importantly, variation in parameter estimation from the 10 data sets is considered in terms of heuristic methods used by experienced electrochemists in different laboratories, and the relevance of the electrode surface state for each of the experiments is examined.
(1)
where A and B are both soluble in the electrolyte-containing solution phase. Even in this case there are usually at least five unknown parameters that need to be specified in the model used for the simulation. These are Ru, Cdl, E0, k0, and α where Ru is the uncompensated resistance, Cdl is the double layer capacitance, E0 is the reversible formal potential, k0 is the heterogeneous electron transfer rate constant, and α is the charge transfer coefficient. This assumes that terms like diffusion coefficient and electrode area are known from independent measurements and that Butler−Volmer theory is employed to describe the electron transfer process. In principle, all the information needed for theory− experiment comparisons is contained in the raw data, and the Fourier Transform (FT) approach represents a mathematical manipulation that is used to aid the experimenter in interpreting the result. The FT method also filters out mains frequency and random noise components. Access to resolved harmonics provides the experimenter with compact bundles of information which have different levels of sensitivity to each of the unknown parameters that are included in the theoretical description. Nevertheless, by definition, the original current− time profile contains all the information (unfiltered) needed to describe the electrode process, so whether an FT or equivalent mathematical transformation is needed can be debated. A major issue is the fact that analytical mathematical solutions for the theoretical model are rarely available, and numerical simulations are invariably needed to solve the model and mimic the experimental behavior. In principle, statistical methods of data analysis are available to assist in theory−experiment comparisons. However, these have rarely been used to date. Rather, empirical heuristic methods are used to decide when a “good” fit of theory and experiment has been achieved. Effectively, this means the experimenter decides on the basis of experience, rather than statistics, when satisfactory experiment−theory agreement has been achieved. In this paper, a comparison of electrode kinetic parameters derived from fully automated and detailed analysis of the raw data and those provided heuristically with the aid of the Fourier transform−inverse Fourier transform modus operandi are considered, along with the role of systematic (mains frequency) and random noise, and the advantages and disadvantages of the variety of data treatment methods available for undertaking theory versus experiment comparisons. In particular, results from 10 consecutive experimental data sets collected for the quasi-reversible reduction of ferricyanide to ferrocyanide
■
EXPERIMENTAL SECTION Materials. Water purified by reverse osmosis (Arium 61315, Sartorius) was used for preparation of the aqueous 3 M KCl electrolyte. Tetra-(n-butyl)ammonium hexafluorophosphate ((n-Bu)4NPF6); 98%, Wako) was recrystallized twice from ethanol (96%, Merck Emplura) and dried under a vacuum at 40 °C prior to use. Ferrocene (Fc; 98%, EGA Chemie), K3[Fe(CN)6] (>99.0%; Sigma-Aldrich), KCl (>99.5%, Merck Emsure), and acetonitrile of Lichrosolv HPLC grade (Merck) were used as received. For deaeration of the electrolytes, high purity nitrogen (99.999%, O2 < 2 ppm) saturated with the vapor of the appropriate solvent was used. Electrochemical Instrumentation and Procedures. Electrochemical experiments were performed in three-electrode cells at ambient temperature (24 ± 1 °C) using an Epsilon Electrochemical Workstation (BAS) or a homemade instrument that has the capacity of being used for Fourier Transformed AC Voltammetric (FTACV) analysis.7 High surface platinum wire separated from the working electrode compartment by a glass frit was used as the auxiliary electrode. The reference electrode was Ag/AgCl (3 M KCl), installed inside a Luggin capillary positioned in the vicinity of the surface of the working electrode and filled with the supporting electrolyte. All potentials are reported versus this Ag/AgCl reference electrode. 11781
dx.doi.org/10.1021/ac4022105 | Anal. Chem. 2013, 85, 11780−11787
Analytical Chemistry
Article
electron transfer (Butler−Volmer theory), mass transport (Fick’s laws of diffusion), Ohm’s law to model Ru, and double layer capacitance theory to describe the background current. Comparison of predictions derived from simulations and experimental data form the basis of estimation of E0, k0, α, Ru, and Cdl in what is known as the solution to the inverse problem. Model. For the purposes of this paper, we will consider a chemically reversible single-electron transfer reaction at a macro-disk electrode as described in eq 1. We proceed as in previous work,28 by considering the reaction
A glassy carbon (GC) disk (ca. 3 mm diameter) supplied by BAS was employed as the working electrode. Prior to each set of measurements, the surface of this electrode was thoroughly polished with alumina powder (Buehler; 0.3 μm) on a wet polishing cloth (BAS). After polishing, the electrode was repeatedly washed with purified water and subjected to sonication (FXP 10M, U-LAB Instruments) in water for 10− 20 s. To provide efficient removal of the alumina powder from the electrode surface, after sonication, the electrode was carefully wiped with a clean wet polishing cloth free of Al2O3 powder, washed with water, and again sonicated for 10−20 s in high purity H2O. Finally, the working electrode was dried in a nitrogen stream. In so far as possible, each electrode preparation was undertaken under identical conditions; however, as will be seen later, the electrode kinetics differ substantially from experiment to experiment. The surface area of the electrode was estimated as 0.070 cm2 from the Randles− Sevcik relationship and data obtained by measuring DC voltammograms for oxidation of Fc at a known concentration in acetonitrile 0.1 M (n-Bu)4NPF6 electrolyte. The value of the diffusion coefficient of Fc in this electrolyte was assumed to be 2.4 × 10−5 cm2 s−1.23,24 Experimental AC voltammetric data for the reduction of 1.0 mM K3[Fe(CN)6] dissolved in an aqueous electrolyte (3 M KCl) were acquired as 10 independent sets of measurements for each solution. For each measurement, the same GC electrode was employed with a freshly polished surface. The aim was to assess the reproducibility of the parameters associated with the [Fe(CN)6]3−/4− process, recognizing that the electrode pretreatment history is not identical and that this is a pathologically variable and notoriously nonreproducible reaction whose response is highly dependent on the state of the carbon electrode surface.10−20 A highly concentrated 3 M KCl electrolyte was chosen for the [Fe(CN)6]3−/4− electrode kinetic studies in order to minimize the IRu (where I is the current and Ru is the uncompensated resistance) drop. To account for this small influence of IRu in simulations, Ru was determined for each experiment by electrochemical impedance spectroscopy at potentials devoid of any Faradaic current using the FTACV instrument. However, it also needs to be noted that the electrode kinetic characteristics of the [Fe(CN)6]3−/4− process, as well as being strongly electrode surface dependent, are also a function of the identity and concentration of the supporting electrolyte cation, as noted by the very comprehensive study by Campbell and Peter25 and references contained therein. The value of the diffusion coefficient of [Fe(CN)6]3−, D, used in simulations of the AC voltammetry was 7.2 × 10−6 cm2 s−1. This parameter was initially determined by analysis of the semi-integrals of background corrected DC cyclic voltammograms, taking into account edge diffusion effects,26,27 and then confirmed by simulations of DC cyclic voltammetric data. Prior to experiments, the glassware was cleaned by soaking in a H2SO4/H2O2 (1:1 vol/vol) mixture and then thoroughly washed with purified water and dried in an oven at ca. 130 °C for 60 min.
kox
A HooI B + e−
(3)
k red
where kox and kred are the heterogeneous rate constants for the oxidation of A and reduction of B, respectively. We assume that convection and migration can be ignored, the former via use of a stationary electrode and unstirred solution and the latter due to an excess of supporting electrolyte, and therefore that mass transport can be modeled by semi-infinite planar diffusion. If we also assume equal diffusion coefficients for A and B, then Fick’s second law gives ∂ci ∂ 2c = D 2i ∂t ∂x
(4)
where for i = A and B, ci is the concentration of A and B, respectively, D is the (common) diffusion coefficient, and x is the distance from a point in the solution phase to the working electrode surface (x = 0). We complete the system with the initial conditions cA(x , 0) =c∞ cB(x , 0) =0
(5)
the boundary conditions as x → ∞ cA → c∞ cB → 0
(6)
and at the electrode surface ⎡ ⎡ (1 − α)F ⎤ ⎡ ∂c ⎤ (E (t ) − E 0 ) ⎥ = k0⎢cA exp⎢ D⎢ A ⎥ ⎣ ∂x ⎦x = 0 ⎣ RT ⎦ ⎣ ⎡ − αF ⎤⎤ (E (t ) − E 0 )⎥⎥ −cB exp⎢ ⎣ RT ⎦⎦
(7)
along with ⎡ ∂c ⎤ If = D⎢ A ⎥ ⎣ ∂x ⎦x = 0 FS
(8) 0
In these equations, k represents the standard rate constant defined at the formal potential E0, α is the charge transfer coefficient, If is the Faradaic current, S is the electrode surface area, F is Faraday’s constant, R is the universal gas constant, and T is the absolute temperature. If the applied potential, Eapp(t), is chosen to be a sinusoidal oscillation of amplitude ΔE and angular frequency ω, where ω = 2πf with f being the frequency in Hertz, superimposed onto a cyclic DC linear ramp, then
■
RESULTS AND DISCUSSION Theory. The implementation of any quantitative form of analysis of the electrode kinetics requires a model that gives rise to predictions of the voltammetric data, for a given set of input parameters. In the case of the mechanism associated with eq 1, the model will be constructed on concepts derived from
Eapp(t ) = min(Estart + νt , Estart + νtmax − νt ) + ΔE sin(ωt ) 11782
(9)
dx.doi.org/10.1021/ac4022105 | Anal. Chem. 2013, 85, 11780−11787
Analytical Chemistry
Article
Table 1. Parameters Recovered from a Simulation with 5% Noise Added, along with the Standard Deviation (std)a experiment 1
a
experiment 2
parameter
simulated
recovered
std
simulated
recovered
std
α E0/V k0/cm s−1 Ru/Ω Cdl/μF cm−2
0.50 0.000 0.1000 2000 25
0.50 0.000 0.1040 2000 25
0.00 0.000 0.0185 15 0.0
0.50 0.200 0.1000 140 500
0.50 0.200 0.1000 140 500
0.00 0.000 0.0051 0.186 0.0
For each parameter, the recovered value reported represents the average from ten runs.
where ν is the scan rate, Estart is the initial potential, and tmax is the end time of the experiment. As a consequence, the effective potential, E(t), is defined as E(t ) = Eapp(t ) − Edrop = Eapp(t ) − ItotR u
current. The approximate current, Inoisy, is consequently given at each time step by Inoisy = Iexact + prImax
(10)
where p represents the percentage of noise, r is a realization from a N(0,1) distribution, Iexact is the numerical solution for the current, and Imax is the maximum current value recorded. We will show later that this is consistent with the noise recovered from experimental data sets. The two experiments referred to above were then repeated, this time including simulated noise at a level of 5%, with the results for the recovered parameters averaged over 10 repeats shown in Table 1. A comparison of parameters used in the simulation and those recovered indicates that even after allowing for a substantial amount of measurement inaccuracy, the algorithm is able to recover the parameters with a high degree of accuracy in a simulation environment, where it is guaranteed that the underlying model is fully representative of the chemistry. Heuristic Approaches. “Heuristic” approaches have been used to process most voltammetric data. In the FT-assisted AC voltammetric method, most commonly an experienced electrochemist will use the experimental data, along with tools designed to simulate the harmonics that would be obtained from a theoretical voltammetric experiment where the five parameters mentioned above are chosen, and follow a recipe to adjust these parameters and attempt to visually match the two. An example heuristic method developed in an earlier study,32 and as outlined in Figures 1 and 2, uses the aperiodic DC component to compute E0 from the average of the reduction and oxidation peak potentials and at least three higher harmonics to match the parameters, followed by a check that all available harmonics satisfactorily match the experimental data. In order to assess potential experimenter bias, we have followed the heuristic recipe independently at both the Monash and Oxford University laboratories to give the results shown in Table 2. As expected in operator-dependent heuristic methods, individuals implemented some variations to the protocols, based on their experience. No consultation took place during the course of these independently undertaken data analysis exercises, so some variation in regard to the parameter values reported is to be expected. Nevertheless, the results show good agreement, with slight variations, particularly in α and in some of the protocols adopted (see footnotes in table). Clearly k0 exhibits significant variation from experiment to experiment. FT Matching. The second FT method we consider incorporates automated harmonic matching, performed by the Nimrod software using parallel processing resources.33 This can allow all parameters to vary and optimizes using a simplex algorithm on the separated harmonics, obtained as done previously.34,35 However, it is still partially heuristic, as a
where Edrop models the loss due to uncompensated resistance, Ru. The total current, Itot, is a combination of the Faradaic current and the background capacitive current, Ic, which may be modeled on the basis of double layer theory, so that Itot = If + Ic = If + Cdl
dE dt
(12)
(11)
with Cdl representing the double layer capacitance. The model is then suitably nondimensionalized,28 and the resulting nonlinear system is solved using a backward Euler finite difference scheme.29 The derivation for a reduction process, such as the [Fe(CN)6]3−/4− process, is analogous, except for the sign of the current which, where appropriate, is negative rather than positive. Global Method of Parameter Recovery. In the fully automated global method, we consider parameter recovery achieved in the time domain, which we will then contrast with outcomes from other methods which are currently used, employing analysis in the frequency domain assisted by Fourier transformation. The optimization in the global approach is formulated as a least-squares problem (see the Supporting Information), and we take as the result the parameters which minimize the difference between simulated and experimental voltammograms under this norm. To solve the optimization problem, we have used a QuasiNewton method30 implemented in the NAG31 C library as e04fcc. This algorithm is designed for unconstrained minimization, but to keep the parameters physical, we have imposed constraints artificially by manually setting the objective function to be orders of magnitude larger should parameters fall outside the given ranges. The simplest test case is to generate simulated data and attempt parameter recovery with no noise added to mimic measurement uncertainties that will be present in experimental data. We have demonstrated this for two examples using extreme values of Ru and Cdl. The parameters used and results showing perfect agreement are provided in Table S1. Other parameters used in the initial exercise were Estart = −0.3 V, Emax = 1.0 V, ΔE = 0.08 V, f = 9 Hz, ν = 0.1 V s−1, D = 2 × 10−5 cm2 s−1, S = 0.071 cm2, c∞ = 1 mM, and T = 297 K. As one would expect, in the absence of noise the algorithm is able to identify the parameters within machine tolerance when an exact match is possible. We then simulated noise and added this to our data. As motivated by ref 8, we have chosen to model this as samples from a Gaussian distribution with zero mean and standard deviation scaled around the maximum 11783
dx.doi.org/10.1021/ac4022105 | Anal. Chem. 2013, 85, 11780−11787
Analytical Chemistry
Article
Figure 1. A recipe for a heuristic method of analysis of AC Fourier transformed voltammetric data as used at Oxford University.
judgment has to be made as to how to deal with ringing and sideband issues caused by the assumption of periodicity when using the Fourier transform, which can be severe in the higher harmonics, as can be seen for example in Figure 2 (fifth harmonic), and also which parameters to fit. Figure S1, showing the fifth harmonic, outlines how this was achieved. In brief, to minimize contributions from ringing, the data for the first and last second of the experiment were ignored. This method is referred to as a computer-assisted heuristic approach. In this form of data analysis, Simplex optimization was performed for the second, third, fourth and fifth harmonics to minimize the objective function LNim, which is calculated as the relative rootmean-square deviation:
Figure 2. The DC component along with the first, third, and fifth harmonics as used in the heuristic approach outlined in Figure 1 and described in the text.
N
L Nim =
∑i = 1 (f exp (xi) − f sim (xi))2
In order to extract the noise from a data set, we make use of a Fourier transform and then look at the power spectrum, where power is defined by
N
∑i = 1 f exp (xi)2
(13)
where fexp(x) and fsim(x) are experimental and simulated functions, respectively, and N is the number of data points. The results from using this method are shown in Table 3. In this case, Cdl was always fixed at the values obtained heuristically at Monash, as the contribution from the double layer capacitance to the relevant harmonics is negligible. Time Domain Matching. The final and fully automated method is time domain matching, as outlined earlier. As this method does not make use of the Fourier transform as a filter, we must include a consideration of the potential impact of noise resulting from measurement error from the instrumentation, which is mostly removed when using the FT-based methods.
Power =
Re(-(Itot))2 + Im(-(Itot))2
(14)
where Re and Im are the real and imaginary components, respectively. Figure 3 shows clearly that the harmonics from the sinusoidal perturbation can be identified, and then these, along with the DC component, can be filtered as “the signal” with everything remaining behind defined as “the noise.” Subsequently performing an inverse Fourier transform returns the noise to the time domain, where it can be analyzed. We have followed this procedure for each of the 10 repeats of the [Fe(CN)6]3−/4− experiment. Figure 4 shows one example of 11784
dx.doi.org/10.1021/ac4022105 | Anal. Chem. 2013, 85, 11780−11787
Analytical Chemistry
Article
Table 2. Parameters Recovered from the 10 repeats of the [Fe(CN)6]3−/4− Experiment Using a Heuristic Approach at Monasha (Mon.) and Oxfordb (Ox.) Universities k0/cm s−1
E0/V
Cdl/μF cm−2
α
exp. no.
Ru/Ω
Mon.
Ox.
Mon.
Ox.
Mon.
Ox.
Mon.
Ox.
1 2 3 4 5 6 7 8 9 10
6 4 5 5 6 7 8 7 8 7
0.214 0.214 0.214 0.214 0.214 0.215 0.215 0.216 0.216 0.217
0.213 0.214 0.213 0.213 0.213 0.213 0.214 0.217 0.217 0.217
0.0105 0.0090 0.0065 0.0055 0.0095 0.0080 0.0095 0.0160 0.0180 0.0115
0.0110 0.0089 0.0067 0.0058 0.0091 0.0080 0.0100 0.0150 0.0190 0.0110
17 18 16 16 18 20 23 24 25 23
17 19 16 17 19 20 24 24 22 23
0.52 0.52 0.52 0.52 0.52 0.52 0.52 0.51 0.51 0.51
0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50
a
A similar protocol to that in Figure 1 was used, but in this case, Ru was taken to be the experimentally measured value and E0 was deduced from peak potentials in odd harmonics and minima in even harmonics. bThe heuristic approach detailed in Figure 1 was used. In this evaluation, Ru was estimated to be
