Integration of Heuristic and Automated Parametrization of Three Unresolved Two-Electron Surface-Confined Polyoxometalate Reduction Processes by AC Voltammetry

The thermodynamic and electrode kinetic parameters that describe each of the three unresolved proton-coupled twoelectron transfer processes of surface-confined Keggin-type phosphomolybdate, [PMo12O40] 3 adsorbed onto glassy carbon electrode in 1.0 M H2SO4 have been elucidated by comparison of experimental and simulated AC voltammmetric data. Modelling of this problem requires the introduction of over 30 parameters, although this may be reduced to about half this number when intelligent forms of data analysis are introduced. Heuristic (i. e., an experimenter based trial and error method) and automated data optimization approaches are integrated in this very extensive parameter estimation exercise. However, obtaining a unique solution remains challenging for reasons that are outlined. In the final analysis and using the automated strategy, estimates of six reversible potentials, lower limits of the six electron transfer rate constants, the double layer capacitance, uncompensated resistance and surface coverage are reported, with others (such as the charge transfer coefficient) present in the model being unobtainable for reasons that are provided. The fit to experimental data using parameters obtained by automated data optimisation is excellent and slightly superior to that obtained by heuristic analysis. The parameters obtained by either method account for differences in shapes and current magnitudes of each of the overall two electron processes.


Introduction
Voltammetric theory for electrode processes comprising an extensive series of coupled heterogeneous electron transfer steps and homogeneous chemical reactions is now very well established. [1][2][3][4] Generation of theoretical data derived from a designated model is known as the forward problem. However, for a complex electrode process, obtaining a large number of unknown parameters that have to be deduced by comparison of experimental and theoretical data, in what is termed the inverse problem, often still remains essentially unmanageable with respect to obtaining both a complete and a unique solution. In essence, solving the inverse problem requires capturing substantial amounts of very high quality experimental data and undertaking repetitive comparisons with theoretical data deduced from a model until acceptable agreement is achieved. When even the simplest possible process in which the oxidized form of the electroactive species (Ox) is converted into its reduced form (Red) as summarized in Eq. 1, and modelling is undertaken assuming that the Butler-Volmer [1,3] relationship applies and mass transport occurs solely by planar diffusion, there is a minimum of 5 parameters that have to be estimated; viz formal reversible potential (E 0 ), heterogeneous charge transfer rate constant (k 0 ), charge transfer coefficient (α), double layer capacitance (C dl ) and uncompensated resistance (R u ), assuming diffusion coefficients (D ox and D red ) and other relevant parameter values are known from independent measurements. [5] If chemical steps are coupled to multi-electron transfer, then on the order of ten unknown parameters will probably need to be quantified. [6][7][8] When addressing a problem of even greater complexity, it may have to be concluded that full parameterization is impossible to achieve when experimental error and model uncertainty are taken into account.
Ox þ e À G E 0 ;k 0 ;a H Red: The forward process of predicting the theory using a proposed model, once also highly demanding when computer coding of all steps was required in each study, can now be achieved routinely with user friendly, commercially available software packages such as DigiSim, DigiElch, KISSA or by using freeware that can be downloaded from the web such as MECSim. [9] Now it is the inverse problem of deciding which model and combination of parameters best describes the experimental data and how good is the fit that usually presents the daunting challenge. Typically, the experimentalist who collected the data may elect to "guess" the model that is applicable and rely on experience to fit the data by essentially trial and error procedures in what almost invariably becomes an extensive series of very tedious theory-experiment comparisons.
In this heuristic approach, the experimentalist decides empirically when an acceptably good fit of data has been achieved and then provides a report of the parameter values that fit the mechanism.
As an alternative to the fully heuristic method, multiparameter fitting aided by computationally efficient data optimization or more sophisticated approaches are now available to assist with solving the inverse problem in voltammetry. [8,[10][11][12][13][14][15][16][17][18] Data optimization methodology, underpinned by statistics and facilitated by high speed computing, has been developed to support many branches of science. Now, just as there are many voltammetric simulation packages available for modelling the forward problem, there is an extensive range of software packages available to support complex theory-experiment data optimization exercises. Nevertheless, in voltammetry, problems requiring an extensive number of parameters to be estimated are rarely attempted. In particular, if the data set available is inadequate as often is the case when using DC cyclic voltammetry (insufficient data points or range of scan rates etc.), then the significance of the outcome of a multi-parameter fitting exercise is likely to remain equivocal. [6,[19][20][21] Furthermore, a range of other difficulties can be encountered, and thus data optimisation methods are far from flawless, and still require heuristic input, as will become apparent in the present study.
Computer aided data optimisation has been applied, at various levels of sophistication and automation in DC voltammetry, to relatively simple problems requiring the estimation of no more than about five parameters (see for example, [8,[13][14][15][16][22][23][24][25][26][27] and references cited therein). Typically in DC voltammetry with automated data optimisation, each experimental data set is interrogated using a particular simulation model and say four parameters such as E 0 , k 0 , α, and D (diffusion coefficient) in the case of a simple solution soluble electron transfer reaction as reported by Scharbert and Speiser. [25] In the minimal and probably most common application, and indeed one available with DigiElch and DigiSim software packages, individual least squares comparisons of differences in experimental and simulated data are undertaken, using only one simulated data set at a time, all being selected by the experimentalist, with lowering of residuals in least squares trends directing the direction of parameter changes. This can be even more tedious than the heuristic method, but offers a scientifically improved level of comparison of experimental and simulated data. In automated data optimisation analysis, the parameter space for all unknowns to be estimated in the model are selected within the expected range and with a selected level of resolution. For example if R u is a parameter of interest, multiple simulations with any value in the range 0 to 100 Ω could be undertaken at 1 Ω intervals as part of the data optimisation exercise. Additionally, a randomly selected initial set of all parameters is computa-tionally selected and the parameter combinations are varied systematically seeking to find the best fit from what can be a very large number of combinations. Consequently, one needs to ensure that computational time is not too excessive by sensibly limiting the parameter range for searching and using Simplex or other well established procedures. These data optimisation approaches are usually referred to as multi-variate or multi-parameter fitting analysis exercises. However, experimenter intervention remains part of the process, as in the heuristic method. With a complex problem of the kind considered in the present study, it is highly risky and not even plausible to simply allow the computer to automatically interrogate every possible combination of parameters without some prior information provided by the experimentalist. Problems with localised minima, over-parameterisation or under-parametrisation may arise and chemically non-acceptable parameter values can then be inadvertently and incorrectly reported. This is where having undertaken a heuristic form of analysis is valuable with the final chemically sensible set of parameters reported being generated by integration of both approaches. All of the above-mentioned limitations and pitfalls will be exposed and discussed in the present study. It of course follows from the above commentary that the tediousness of the use of the heuristic methods should be kept to a minimum by use of judicious strategies and that the parameters ultimately deduced from the automated data optimisation procedure should be relevant (sufficiently sensitive) to the simulation model and be superior to those selected from the purely heuristic approach.
In recent work, the Oxford and Monash University Groups have been developing protocols to address the issues that arise when attempting to parametrize increasingly complex mechanisms. In summary, very large data sets containing extractable variable time (frequency) domain information are now collected at high resolution using instrumentation having 18 bit DAC and ADC converters. [28,29] The waveforms developed for these studies are based on employing a large amplitude periodic waveform superimposed onto a DC potential ramp. The use of the Fourier transform algorithm as part of the protocol associated with Fourier Transform AC voltammetry allows a series of harmonics to be resolved as well as the aperiodic DC component. Thus, features related to electrochemical impedance spectroscopy and DC cyclic voltammetry are made available simultaneously as well as additional ones. [26,[28][29][30][31] Experiment-theory comparisons have then been undertaken at levels ranging from fully heuristic to multi-parameter fitting in attempts to uniquely define the thermodynamic, kinetic, mass transport, capacitance and resistance related parameters that are included in the model used to generate the simulated data. This modelling approach uses parameters that have a direct relationship to the physical chemistry associated with the reaction mechanisms, unlike the use of equivalent circuit models traditionally employed [3] in electrochemical impedance spectroscopy (EIS), although both approaches are of course mathematically interchangeable. [16,32,33] In the present study, we have taken advantage of access to ever expanding computing power as well as software available from many sources to quantify a complex mechanism which even in a simplified form contains 17 unknown parameters in the model. These parameters were quantified initially by both fully heuristic and fully automated data optimization methods. However, ultimately, a hybrid approach in which an automated data optimization strategy is informed by knowledge gained from the heuristic method of data analysis was found to provide the best statistical fit of theoretical to experimental data.
Polyoxometalates (POMs) which are of interest in this paper have been widely employed in chemistry in diverse fields such as electrocatalysis and photocatalysis, [34,35] sensors, [36,37] batteries and capacitors [38,39] and are of interest in many branches of science and technology. Since many practical devices based on POMs exploit their extensive redox capabilities, detailed studies of their electrochemistry are needed to facilitate their development. The electrochemistry quantified in this study is the six electron reduction of the surface confined polyoxometalate [PMo 12 O 40 ] 3À . This inorganic cluster exhibits a Keggin type structure ( Figure 1a) which contains 12 molybdenum atoms in oxidation state VI that can be reduced to mixed valent forms in multi-electron steps to give highly negatively charged and very basic mixed-valent products that facilitates the coupling of electron and proton transfer reactions.
The DC voltammetry of the [PMo 12 O 40 ] 3À anion, at solid electrodes and polarography at the dropping mercury electrode have been extensively reported when dissolved in aqueous electrolyte media, molecular organic solvents containing supporting electrolytes or ionic liquids. In aqueous media, [PMo 12 O 40 ] 3À is known to spontaneously adsorb onto glassy carbon, [40] gold, [40] silver, [40] boron doped diamond [41] and reduced graphene oxide electrode surfaces. [42] DC and AC voltammograms obtained at a glassy carbon (GC) electrode in 1.0 M H 2 SO 4 are displayed in Figures 1b and 1c respectively. The three very well-defined surface confined reduction processes evident in these voltammograms are designated as Processes I, II and III. Each represent overall two electron-two proton coupled processes that can be represented by Eq. 2-4. [43] ðProcess IIÞ ð3Þ The peak potentials are located at about 350 (Process I), 225 (Process II), and 10 (Process III) mV vs Ag/AgCl in 1.0 M H 2 SO 4 , but depend on acid concentration, as expected if a net twoproton transfer reaction accompanies a net two-electron transfer reaction. Even more extensive reduction occurs at more negative potentials than displayed in Figure 1 to give an analogous Process IV (not shown) but this leads to rapid dissolution of solid on the voltammetric timescale. Use of lower acidities also facilitates dissolution of reduced forms of the POM. Furthermore, the designation of the fully oxidized [PMo 12 O 40 ] 3À as completely unprotonated in 1.0 M H 2 SO 4 is unlikely to be correct. Further points for noting which can be gained from perusal of Figure 1 and which need to be accommodated in simulationexperiment comparisons are that the first two processes partially overlap, Process II has a larger peak current magnitude than either Processes I or III and the shapes of all three process differ. Importantly, for each of the three reduction steps it has been proposed [4,41,[44][45][46][47] that two unresolved reversible one electron transfer steps of purely adsorbed material occur as in Eq. 5-6 where Additionally, C dl , R u , and the adsorption isotherm need to be modelled, as in principle do the thermodynamics and kinetics of the chemical (acid-base) reactions coupled to electron transfer. In studies of solution soluble redox couples, control experiments can be undertaken at the bare electrode which allows C dl and R u to be measured independently of the faradaic process. However, chemical modification of the electrode alters the background current, meaning that C dl and R u values obtained at the bare electrode are not available for use as independently known parameters. Thus, in principle, there are in excess of thirty parameters present in the full POM reduction modelling exercise. As indicated above, even at this level of complexity, all parameters could be included in the forward problem using MECSim, DigiElch or other software packages. However, it is not likely to be realistic to simply undertake a greater than 30 parameter determination exercise using a manual heuristic approach, or even an automated approach to the optimization experiment-theory fitting exercise (i. e., solving the inverse problem), and expect to obtain a chemically and statistically unique value for each one. In a practical sense it will be shown in this study that the best opportunity to approach this goal is by integration of knowledge gained from both heuristic and automated data optimization approaches to achieve a high level of agreement between experiment and theory by parameterization when there are many parameters that could significantly influence the voltammetry.
Finally, in this Introduction it is emphasized that the data obtained experimentally has to be sufficient in quality and quantity to justify the conclusions reached in a parameterization exercise. AC voltammetry in square wave, sinusoidal or other forms as well as EIS provide significant advantages over traditional DC methods in quantitative studies of electrode processes. [29,31] If large rather than traditional small amplitude sinusoidal AC signals are employed, the significantly amplified higher-order AC harmonic components become available that are virtually free of background capacitance current and are also highly sensitive to electrode kinetics. The ability to resolve and individually analyse the aperiodic DC, fundamental and higher order harmonics can be crucial in heuristic forms of data analysis. For simple problems, comparison of experimental data with the numerical simulations of appropriate models have provided good estimates of E 0 , k 0 , α, R u and C dl from a single experiment, using both heuristic and data optimization forms of experiment-theory comparison. Large amplitude Fourier transformed AC voltammetry (FTACV) has been applied to the determination of the electrode kinetics of surface-confined enzymes. [48][49][50][51] In some cases, enzyme and metalloprotein voltammetry mimics that found with [PMo 12 O 40 ] 3À in the sense that acid-base chemistry is coupled with an unresolved 2 electron transfer process to provide a mechanism and hence data analyis problem akin to a combination of Eq. 2-4, 5 and 6. In this study, experience acquired in parameterization of inherently simpler systems is crucial in addressing the vastly more complex problem of parameterization of the initial six electron reduction processes associated with the surface confined [PMo 12 O 40 ] 3À in 1.0 M H 2 SO 4 at a GC electrode.

Chemicals and Reagents
Sulfuric acid (H 2 SO 4 , Merck,) and phosphomolybdic acid (H 3 PMo 12 O 40 , Sigma-Aldrich, � 99.99 %), were used as received. Water employed for preparation of all aqueous solutions was obtained from a Milli-Q water purification system.

Instrumentation and Procedures
AC voltammetric measurements were undertaken using in-house instrumentation [28] with a sine wave of amplitude (ΔE) of 20 mV half-peak-to-peak and frequency (f) of 9.02 or 60.05 Hz being superimposed onto the DC ramped potential which had a known scan rate (v). The amplitude of 20 mV is smaller than used in most of our previous studies because we wish to minimise overlap of processes I and II. This means we are restricted to detection of only four AC harmonics with good signal to noise ratios. DC voltammetric data also were obtained with this instrument. A conventional three-electrode electrochemical cell under potentiostatic control was used to collect experimental data at 25°C. GC (diameter = 3.0 mm, CH Instruments, USA) was used as the working electrode, Ag/AgCl (3 M KCl) as the reference electrode and platinum wire as the auxiliary electrode with 1.0 M H 2 SO employed as the supporting electrolyte. To prepare POM modified electrodes, the GC surface was polished with a 0.3 μm alumina aqueous slurry on a polishing cloth, rinsed with water, sonicated and rinsed again with water. The cleaned GC electrodes were then dried under a stream of nitrogen gas and placed in 0.1 mM phosphomolybdic acid (H 3 PMo 12 O 40 ) aqueous solution for 1 min. The resultant POM modified electrodes were rinsed several times with water to remove unbound or loosely bound solids before being used in electrochemical experiments.  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55 56 57

Simulations and Data Analysis
Simulations of total current AC voltammograms used in initial heuristic forms of data analysis were carried out with MECSim (Monash Electrochemistry Simulator) software. [9] In the automated data optimisation studies, modelling was undertaken using the PINTS (Probabilistic Inference on Noisy Time Series) software [52] and code specifically written for use with this method. Details of the model and parametrization employed with this approach are provided in the Results and Discussion section. In the FTACV version of data analysis, the experimental and simulated total current time domain data were converted to the frequency domain to generate the power spectrum. Frequencies corresponding to the AC harmonics and the aperiodic DC component were then selected from the power spectrum and this data subjected to band filtering with a rectangular window. Finally, inverse Fourier transformation was employed to obtain the resolved DC and AC components as a function of time. An estimate of the uncompensated resistance R u used in the heuristic form of data analysis was obtained experimentally from the R u C dl time constant [3] at potentials where no faradaic current was present. In total current AC voltammograms, the full data set are provided in the figures. However, FTACV data are presented in the envelope form to facilitate visual comparisons of AC harmonic based experimental and theoretical results

The Model
To simulate the total current AC voltammetric data, a mathematical model has to be employed. To mimic the 3 surface confined processes summarized in Eq. 2-4, a series of well-known relationships were employed in combination to give a total model as follows: 1. Electron transfer model: the Butler-Volmer relationship was used for each of the six electron transfer steps which requires the introduction of six E 0 , six k 0 and six α parameters. Alternatively, the mathematically more sophisticated Marcus-Hush relationship could have been used to model the electron transfer steps and thermodynamic and kinetic dispersion could have been included to reflect any heterogeneity in electron transfer associated with nonuniformly surface confined POM. 2. Acid-base chemistry: proton transfer reactions coupled to electron transfer were assumed to be diffusion controlled and hence reversible and subject to a solely thermodynamic description. On this basis, the unknown acid (pK a ) equilibrium constants have been combined with the E 0 values so effectively only six reversible potential thermodynamically relevant parameters with their six charge transfer rate constants and charge transfer coefficient values need to be parameterized. This approach means that the true individual E 0 , k 0 , α, and pK a values remain unknown. However, for convenience, and as is standard practice in most studies, the usual E 0 , k 0 and α notation is retained in presentation of results below. 3. Adsorption isotherm: the Langmuir model was used, which means that a surface coverage parameter is needed. Other isotherms are available such as Frumkin or Temkin. Use of the Langmuir isotherm implies that the POM remains fully surface-confined for the duration of the experiment. That is no dissolution occurs, allowing its diffusion to be neglected in the modelling. 4. Uncompensated resistance: Ohm's law used (R u parameter introduced). 5. Background current: Modelling undertaken assumes that a simple R u C dl time constant applies at each potential and that C dl is independent of potential. In fact C dl depends on potential, particularly in the DC component (Figure 1b), with a lower dependence being evident in the total current ( Figure 1c) or fundamental harmonic (see discussion below) when the DC term is removed. Close to ideal behaviour is evident in the second, third and fourth AC harmonic components, which as predicted theoretically are entirely devoid of background current (see discussion below).
Presumably, the background current is not fully capacitive, particularly on the DC time scale (pseudo-capacitive behaviour). The background current at a POM sub-monolayer covered GC electrode is derived from GC containing functional groups such as quinones etc and POM species whose comparison depends on potential. Thus the POM modified electrode is highly heterogeneous and hence difficult to fully model. 6. Many other parameters are present in the model (see theory below), but a significant number such as AC amplitude and frequency, DC scan rate, electrode area, temperature etc. are assumed to be accurately known.
In summary, the model and protocols chosen for simulation of AC voltammograms left a total of 17 parameters to be estimated by automated data optimization, but less heuristically for reasons given below, under circumstances where there are assumptions introduced along with probable imperfections in several aspects of the model, as noted above.
The sequence of six electron transfer steps (Eq. (5) and (6) where the forward k red and backwards k ox reaction rates are given by the Butler-Volmer relationships The applied potential E r (t) is the input value E(t) minus the Ohmic IR u potential drop associated with the uncompensated resistance R u .
and the input potential is the sum of DC and AC components where ν = 104.3 mV s À for experiements at 60.05 Hz, E s = 600 mV, and E r = À 100 mV.
The ordinary differential equation (ODE) governing the behaviour of θ is given by where the stoichiometry matrix K is given by and the vector c arises due to the 1 in the proportion of G, and is given by We use a backwards Euler discretisation of the time gradient in Eq. 19 with the exception that the total current was treated explicitly. That is, I n tot was used to calculate K. This leads to the following linear system to be solved at each timestep where δt is the timestep, and θ n is the proportion vector at timestep n.
The total current measured is the sum of the capacitive I c and faradaic I f components which are given by where the change of charge e with time is given by The theory for a simpler case involving only two unresolved transfer steps has been considered in the context of FTACV [48,49] as well as other AC methods. [41] Details of methods employed in simulations in voltammetry are available in references [2,3,53]

Heuristic Method of Parameter Estimation
Implementation of the heuristic or exclusively experimenter based trial and error method in a tractable form requires simplifying the problem and obtaining initial estimates of parameters in a sensible stepwise fashion by relying on experience gained with simpler problems. In this exercise, the solution to the forward problem was obtained from the simulation package MECSim with an initial guess made for each parameter. Simulations with new sets of parameters were then iteratively changed in the direction dictated by the experience of the experimenter until this individual decided that an acceptably "good fit" to data had been achieved. Of course a different experimentalist may decide that the "good fit" is achieved with a different set of parameters!
The initial scrutiny of experimental data involved visual interrogation of AC voltammograms obtained from 20 mV amplitude data sets at 9.02 and 60.05 Hz in total current (see Figure 1c for 60.05 Hz), power spectra and the first 4 AC harmonic formats (see Figures in experiment versus experiment comparisons presented below) and a DC cyclic voltammogram (Figure 1b). This overview allowed initial guesses or independent estimates to be made for some parameters along with some deductions as to the validity of the model employed as follows: (a) It was noted in the fundamental AC harmonic that a potential (time) region exists prior to onset of faradaic current for process I that is close to potential independent and was assumed to be purely capacitive. This region was analysed in terms of the R u C dl time constant to give an estimate of R u as 49 + /À 5 Ω. (b) By matching experiment and theory for the same potential (time) region in the fundamental harmonic that is devoid of faradaic current as described in (a) and effectively use of Eq. 25, C dl was estimated to be 7 + /À 1 μF cm À 2 . However, the sloping baseline in the DC cyclic voltammogram and mismatch of positive and negative potential regions of fundamental harmonic devoid of faradaic current, and indeed the aperiodic DC component of FTACV data, it is clear that C dl is potential dependent, rather than potential independent, as assumed in the model. Thus, treating C dl as a potential independent parameter in an experiment versus theory exercise as undertaken in this study will represent an imperfect approximation. More complex potential dependent estimates of C dl could have been included in the model (Eq. 25). However, perusal of third and fourth AC harmonics reveals no background current in these responses, allowing the electrode kinetics to be estimated from these data that are devoid if capacitance current. Consequently, use of the more complex capacitance model was not regarded as necessary, particularly under conditions where the Ohmic drop or I total R u is small, as applies in this study. (c) The surface coverage was found to be G = 30 + /À 5 pmole cm À 2 based on integration of DC current time-data to give the charge associated with reduction of the surface confined POM and use of Faraday's Law. Since individual fully resolved one-electron charge transfer processes are unavailable, the charge used in the estimation of G was calculated from the sum of that derived from processes I and II (n = 4) or process III (n = 2) after background (capacitance current) correction of a DC linear sweep voltammogram. (d) Since all three processes I, II and III are known to involve the overall transfer of two electrons, are sharp and give a series of well-defined higher order AC components, the k 0 values were assumed to be large. The initial guess was that all k 0 values are sufficiently large that with a low frequency of 9.02 Hz that all one electron processes could initially be assumed to be reversible. This is consistent with the square wave studies of Molina et al. [41] who assumed complete reversibility applies at boron doped diamond electrodes with frequencies in the range of 5 to 10 Hz. Furthermore, the shapes and other attributes of each of unresolved one electron transfer steps in processes I, II and III clearly required that the second electron transfer step occurred at a similar (in the case of process II) or even more negative reversible potential than that of the first electron transfer step. To obtain more detail, Figure 2 was generated to provide simulated fundamental and third harmonic AC data where the separation of the two E 0 values in process II fixed at zero (largest current magnitude) while processes I and III had variable separations in the E 0 (ΔE 0 ) values from 0 to 55 mV with all k 0 and α values being set at 1 × 10 8 s À 1 and 0.50 respectively to mimic a reversible process and other parameters set at estimates approximating those deduced in (a), (b) and (c) above. Potentials greater than 55 mV were  Figure, where 01-00-01 means that ΔE 0 = 10 mV for process I, 0 mV for process II and 10 mV for process III. deemed unreasonable. By eye, and with some fine tuning, ΔE 0 = 30, 0 and 27 mV for processes I, II and III respectively gave a satisfactory match of theory and experiment for the 9 Hz low frequency data set. To confirm that the assumption of reversibility is valid and what this means, Figure 3 was generated with simulated fundamental and third harmonic response for the pairs of k 0 values always the same but ranging from 1.0 × 10 À 3 to 10 8 s À 1 with α being set at 0.50 and ΔE 0 = 33, 0 and 45. At the low frequency of 9.02 Hz, all three processes are insensitive to k 0 electrode kinetics in excess of 10 2 s À 1 . Thus it was concluded from this trial and error analysis that ΔE 0 = 30, 0 and 27 mV were appropriate for processes I, II and III which translates to E 0 values of 368, 338, 227, 227 11 and À 16 mV Ag/AgCl and that some, if not all, k 0 values are larger than about 10 2 s À 1 . (e) Examination of the higher frequency 60.05 Hz AC data set implied that the electrode kinetics could possibly be determined from theory-experiment comparisons of their higher order AC harmonics However, the only reasonable trial and error heuristic approach available was to set the k 0 value for each pair of processes to be equal to each other, set α = 0.50, and assume E 0 values and other parameters are known from 9.02 Hz data. Now only pairs of k 0 values and G remain as unknowns. After much tedious trial and error examination of many simulations with various k 0 values along with fine tuning of surface coverage and uncompensated resistance values, the excellent match of theory and experimental data shown in Figure 4 was achieved for AC harmonics 1 to 4 derived from the FTACV protocol with use of simulation parameters contained in the captions to this Figure. Nevertheless, the k 0 values employed in the simulation in Figure 4 should not be regarded as true values. The heuristic analysis essentially only implies that very fast electrode kinetics approaching the reversible limit is evident for all six electron transfer steps included in processes I, II and III as any k 0 value  greater than about 5 × 10 2 s À 1 generates essentially the same simulated data set as when using the parameter values given in caption to Figure 4. In order to be a tractable exercise, the fit to experimental data obtained heuristically had to be achieved with a number of assumptions that limited the number of parameters needing to be estimated. Thus, for example it had to assumed that three pairs of k 0 values are equal. It is now of interest to establish what outcome could be achieved via multi-parameter data optimization and how this exercise should be implemented with 17 unknown parameters and without the assumption of equal k 0 values and independently known surface coverage and uncompensated resistance values. It will emerge that the automated data optimisation approach is still not fully straightforward to implement and that in practice a hybrid approach based on using knowledge obtained above heuristically, supported by computationally efficient data optimization strategies are needed to obtain chemically sensible outcomes.

Parameter Estimation using Automated Data Optimization Methods
The set of parameters to be found in the automated data optimisation is given by the vector P Heuristic data analysis revealed that all electron transfer processes are reversible or very close to reversible and hence very insensitive to the value of α. This is confirmed by automated parameter estimation when α 1 and α 2 are included in the list of parameters to be estimated. However, in order to save computational time, values of α 3 , α 4 , α 5 and α 6 were set at 0.50 rather than estimated by automated data optimisation, although in practice all α values are in fact concluded to be non-determinable.
In the automated data optimisation approach, we also introduced an experimental current trace given by a vector I exp , where each element of the vector I exp i is the measured current at sample time t i . Using the simulated model given previously, we can calculate an equivalent simulated current trace I sim (p), with elements I sim i calculated at the same set of sample times t i . Since the timestep δt is unlikely to match the experimental sample times, we use linear interpolation to resample the simulated current trace to the experimental sample times.
A traditional optimisation function for calculating the optimal parameters given I exp is the 2-norm distance metric FðpÞ ¼ jjI exp À I sim ðpÞjj ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 N Minimising this distance will give an optimal set of parameters. We perform the minimisation using CMA-ES, [54] an evolutionary algorithm for non-linear optimisation, using the following parameter bounds: where R h u and C h dl and G h are the values of R u , C dl and G obtained by the heuristic method.
For this problem, regularisation is necessary to ensure convergence to the optimal parameter set. We find in practice that using the standard 2-norm distance metric often results in the algorithm stalling in parameter regions where the values of the reversible potential E 0 i for the six different reactions have swapped their order. That is, for example, E 0 1 > E 0 3 , even though we know that Process I occurs before Process II.
We regularise the problem by using a minimisation function inspired by the statistical Bayesian framework. Assuming independent Gaussian measurement noise at every timestep, we can write down the log-likelihood of I exp occurring given a parameter set P as [55] LðpÞ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 2s 2 where σ is an unknown measurement noise level that we add to the parameter vector P.
We encode our knowledge of the parameters in a Bayesian normal prior for each of the E 0 i parameters with a mean m i and standard deviation σ i .
The mean encodes our "best guess" for the value of E 0 i and the standard deviation encodes how confident we are in this "best guess". For the remaining parameters we use a uniform prior between the lower and upper bounds given in Eq. 29. Since these uniform priors are independent of the parameter vector P they do not contribute to the minimisation.
Given the large number of experimental data samples provided (N = 32728), there is a danger that the log-likelihood function will dominate over the prior (i. e. the prior will have no effect on the performance of the algorithm), therefore we also scale the log-likelihood based on N. Combining the prior with the scaled log-likelihood gives the proposed objective function to be minimised We use CMA-ES to minimise Eq. 32, with ðm 1 ; m 2 ; m 3 ; m 4 ; m 5 ; m 6 Þ ¼ ð353; 353; 227; 227; À 2:5; À 2:5Þ mV ð33Þ and σ 1 = σ 0 , where σ 0 is allowed to vary between 1 30 ðE s À E r Þ � s 0 � 21 30 ðE s À E r Þ. We therefore have a single "tuning knob" for our automated fitting algorithm, σ 0 , which encodes how confident we are in our guesses for the E 0 values. Figure 5 summarises the results for the automated fitting while varying the width σ 0 of the E 0 priors, using the quasireversible electron transfer model described earlier. For each value of σ 0 , 20 separate runs of the CMA-ES algorithm were done, each resulting in a minimum F score describing the quality of the fit (a lower score correlates in a better fit to the experimental data), along with a corresponding set of optimal parameters P. In total, 190 independent fits needed to be performed, and we used the University of Oxford's Arcus Phase B cluster (comprising of Dual Haswell CPU nodes which have 16 cores per node and a minimum of 64Gb of memory) to perform this analysis over a 12 hour period.
For each value of σ 0 , the final minimum F score for each of the 20 fits is shown in Figure 5 as blue dots (note that many of 20 runs obtain a similar score and thus many of the dots overlap). The number of "good fits" for each value of σ 0 was then calculated, where a good fit is determined as one whose F score was within 1 % of the lowest F for each σ 0 . The number of "good fits" for each σ 0 is shown in the Figure as red crosses, and were only calculated for values of σ 0 where at least one fit was determined (visually) to represent a satisfactory fit to the experimental data. It was found that for σ 0 > 250 mV, none of the 20 optimization attempts resulted in a good fit and therefore the number of "good fits" is not shown for this range.
The results show that a σ 0 < 100 mV is required for an 100 % reliability for the automated fitting (i. e., with all 20 fits resulting in the same satisfactory fit to the experimental data). These results are very promising, since this indicates that an initial guess on each E 0 value only needs to be accurate within 14 % of the potential sweep range for automated fitting to be reliable. An experienced electrochemist can easily read this estimate off a simple current-potential plot of the experimental trace. Figure 6 show a comparison of the best fit with σ 0 = 100 mV, against the experimental data. The first four harmonics are shown for comparison, and it can be seen visually that the fit due to the automated process is excellent. Noteworthy is the fact that best fit k 0 rate constants (see Figure 6 caption) identified for each of the six electron transfer steps present in the model are all very large (greater than 4000 s À 1 ), but the charge transfer coefficient α could not be determined for any process. These features are consistent with the electrode kinetics being extremely close to the reversible limit and probably indistinguishable from reversible processes within experimental uncertainty, as concluded from use of the heuristic method. However, the reversible potentials deduced from the automated data fitting exercises are considered to be superior to those estimated heuristically and reveal that a crossover in their order accounts for the larger peak current and different shape of the second overall two electron process II, although differences in individual reversible potential from the two steps that contribute to each overall two electron processes I,II and III are always small (tens of mV magnitude). In summary, we have shown that an automated fitting process can be reliable, given relatively weak prior information by the user, and result in a fit that is superior to the overly laborious, heuristic method (see evidence provided below).
The large k 0 values found by both the heuristic and automated fitting process indicate that all the reactions are Figure 5. Results of automated fitting of model to experimental data based on the quasi-reversible electron transfer model described in the paper. The left vertical axis shows the objective function F in Eq. 32 versus the scaled standard deviation σ 0 of the normal prior on the reversible potentials (i. e., the range over which the user expects the optimal E 0 i to be found). A higher σ 0 indicates a wider range, or less confidence in the initial guess. The right vertical axis shows the number of good fits out of 20 attempts, where a good fit is defined as one that results in an minimal F within 1 % of the best of the 20 attempts (note, the exact equivalence of estimated parameters in some runs, within the resolution of this plot, accounts for the lower number of data entries that would otherwise be anticipated). The results show that σ 0 � 100 mV is required for automated fitting with a reliability of 100 %. operating close to the reversible limit. Therefore it is instructive to consider fitting an easier model with each k 0 i parameters set to a constant value of 10 4 s À 1 (increasing k 0 i above this value does not change the results significantly since the simulated model is in the reversible limit). Similar results comparing the number of good fits with varying σ 0 are shown in Figure 7. Here we see that the prior information necessary for the reversible potential parameters is much the same as for the quasireversible case, even though we are fitting a far lower number of parameters (9 instead of 17). This indicates that it is the nonlinear response of the system in response to changes in E 0 i that causes difficulty in the minimisation process, rather than the higher dimensional parameter space.
We also provide the best fit for the reversible case in Figure 8. This is, as expected, virtually indistinguishable from the quasi-reversible case, indicating that the reactions are proceeding close to the reversible limit.
In concluding the report on automated data optimisation it is useful to summarise how the substantial number of issues identified in the example studied in this paper were resolved. 17 unknown parameters were present in the model finally selected for interrogation by automated data optimisation with the ranges and resolution of the parameter space searched for each one summarised in Eq. 29. Initially, little constraint was placed on the order or values of any the reversible potentials except that they had to lie between the initial and final values use experimentally. However, chemically impossible values were generated as outputs from this form of computational approach. Only when knowledge gleaned from heuristic analysis that each set of reversible potentials have to be similar and that the logically predicted order set by processes I, II and III has to be followed were sensible parameter estimates obtained. The reversible potentials presented were by far the most difficult parameters to estimate by automated data optimisation. R u for example could be straightforwardly defined to lie in the range 0 to 500 Ω knowing that a value of 49 Ω had  . The left vertical axis shows the objective function Eq. 32 versus the standard deviation σ 0 of the normal prior on the reversible potentials (i. e., the range over which the user expects the optimal E 0 i to be found). A higher σ 0 indicates a wider range, or less confidence in the initial guess. The right vertical axis shows the number of good fits out of 20 attempts, where a good fit is defined as one that results in an minimal F within 1 % of the best of the 20 attempts (note, the exact equivalence of estimated parameters in some runs, within the resolution of this plot, accounts for the lower number of data entries that would otherwise be anticipated). The results show that σ 0 � 100 mV is required for automated fitting with a reliability of 100 %. been determined independently and confirmed at about 50 Ω heuristically. Indeed the finally value reported value of 40.3 Ω is chemically sensible, noting that the AC response is not strongly dependent on the value of this parameter. The surface coverage details also could be set within a range guided by the coulombically determined and heuristically supported values. However, it should be noted that the R u and surface coverage values determined by independent methods were undertaken on different experimental data to those analysed by data optimisation, and need not exactly represent the values relevant to the AC data set. The experimentalist also had to assess the validity of the model used to simulate experimental data. Heuristic analysis work suggested that either quasireversible or fully reversible (implies over parameterisation) models of electron transfer could have been appropriate. Restrictions on α lying in range 0.4 to 0.6 or else set at 0.5 were also made. Due to lack of sensitivity issues it eventuated that only a lower limit on k 0 values could be placed and α could not be estimated as the AC data are very insensitive to their value. In essence the reversible model defines the electrochemistry to a very good approximation. Other modelling issues that had to be resolved were as to whether to use a constant potential or potential dependent value of C dl with the former model possibly representing under-parameterisation. The simpler model, on basis of heuristic examination of data was regarded as acceptable. With respect to the AC data set, data optimisation could have been undertaken on the total current or on the individual harmonics after filtering out say the DC and fundamental harmonic components which contain by the largest contribution from imperfectly modelled charging current. Since the heuristic method taught us that the potential independent C dl assumption should be reasonable as the charging current is not a major contributor to the total current, which is mainly faradaic in nature, we chose to use the total current in the automated data optimisation exercise instead of modelling individual AC harmonics, but carefully checked that each harmonic provided an excellent fit to experimental data. The above choices were guided by prior knowledge held by the experimentalist after undertaking the heuristic form of data analysis, but also took into account computational time which is likely to be excessive with a multi-parameter evaluation exercise of the present kind. Indeed, with respect to the latter issue, a successful parameter evaluation for the reduction of the polyoxometalate would probably not have been possible with the computational facilities available to support this study without the combination of placing significant restrictions on the parameter space searched (equation 29) and having a strategically planned computationally efficient modus operandi.

Conclusions
The voltammetry of [PMo 12 O 40 ] 3À adsorbed onto a GC electrode in contact with 1.0 M H 2 SO 4 consists of three overall two electron-two proton coupled reactions over the potential range employed in this study. Simulations of the six one electron transfer reactions and coupled proton steps associated with the voltammetry is relatively straightforward. However, even with a judicious selection of assumptions that remove the need to include the acid-base chemistry, the simulations requires the input of more than 17 parameters whose values are unknown.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57 estimated, each having a k 0 value in excess of about 2000 s À 1 . In addition, the electrode surface coverage and uncompensated resistance values have been estimated as part of the data optimisation exercise. The values of the two E 0 values associated with each pair of electron transfer reactions that contribute to the overall processes I, II and III are similar which accounts for variation shapes and peak current magnitudes in DC and AC total current voltammograms (see Figure 1). As might be expected, the agreement between theory and experiment is superior when using the automated data optimisation strategy. A convenient way of comparing the level of agreement of experiment and theory achieved with the different data analysis methods is to overlay the experimental AC total current versions and experimental predictions as has been done in Figure 9. Agreement of simulated and experimental data is almost perfect in the data optimisation approach, but is still very good for the heuristic method. Figure 9 also contains a comparison of the theoretically predicted response obtained heuristically and by data optimisation and they are indeed very similar. Clearly the automated data optimisation method is superior with respect to parameter estimation and the outcome can be reproduced in all laboratories. However, at present and unlike when addressing the forward problem where MECSim and other software packages are available, web applications that are specifically designed for addressing the inverse problem in electrochemistry at a sophisticated and automated level are not yet available. In the future we plan to develop such a web application for voltammetry, based on one available for cardiac electrophysiology. [56] Data and Software Availability The experimental data and software used for the automated data optimisation is available at https://github.com/martinjrobins/pom_project.