Anal. Chem. 2006, 78, 1807-1816
High-Dimensional Model Representation of Cyclic Voltammograms Lesław K. Bieniasz*,† and Herschel Rabitz*
Department of Chemistry, Princeton University, Princeton, New Jersey 08544
Digital simulation costs present an obstacle on the way to high-speed, real-time, on-line theoretical analysis of experimental data in cyclic voltammetry. To overcome this difficulty, we propose to use solution mapping based on a correlated, hierarchical expansion of multivariate functions, known as high-dimensional model representation (HDMR). The nonlinear dependencies of the simulated voltammograms on multiple model parameters are represented in the form of compact look-up tables, from which approximate voltammograms can be calculated rapidly by interpolation, for any model parameter combinations from a predefined domain. Most importantly, the HDMR does not suffer from the problem of the exponential growth of the look-up tables with the number of model parameters. The creation of a solution map requires a single effort of simulating many voltammograms. However, once the map is prepared, it can be stored and reused many times without the need to repeat costly simulations. HDMR maps are created and examined for five examples of cyclic voltammetry models at planar macroelectrodes in a one-dimensional spatial geometry under pure diffusion transport conditions. The usefulness of the maps for rapid visualization and exploration of the effects of the parameters on the voltammograms and for rapid simultaneous estimation of many parameters from cyclic voltammetric data is demonstrated through computational experiments. Cyclic voltammetry (CV) is one of the most powerful and commonly employed electrochemical transient techniques, used both for analytical purposes and for kinetic/mechanistic studies in electrochemistry.1 However, the versatility of CV comes at a relatively large cost (compared to some other transient methods) of obtaining theoretical solutions describing CV experiments. Lowcost computable analytical solutions of the relevant initial or initialboundary value problems do not normally exist, except for a few simple kinetic cases. Therefore, computer-aided modeling2 by * Corresponding authors. E-mail:
[email protected]; nbbienia@ cyf-kr.edu.pl, http://www.cyf-kr.edu.pl/∼nbbienia. † On leave from: Department of Electrochemical Oxidation of Gaseous Fuels, Institute of Physical Chemistry of the Polish Academy of Sciences, ul. Zagrody 13, 30-318 Cracow, Poland. Tel. /fax.: (+48 12) 266 03 41. (1) See, e.g.: Speiser B. In Encyclopedia of Electrochemistry; Bard A. J., Stratmann M., Unwin P. R., Eds.; Wiley: New York, 2003; Vol. 3, Chapter 2.1, p 81. (2) See, e.g.: Mocak J. In Encyclopedia of Analytical Science; Worsfold P. J., Townshend A., Poole C. F., Eds., Elsevier: Oxford, 2005; p 208. 10.1021/ac051373r CCC: $33.50 Published on Web 02/04/2006
© 2006 American Chemical Society
digital simulation methods3 is necessary. Despite the considerable progress in the simulation methodology achieved in recent decades, and the systematic increase in computer speed, the simulation of cyclic voltammograms still requires significant time compared to the currently achievable experimental data acquisition times. Particularly long simulation times occur in the case of multispecies electrochemical reaction systems requiring one- or higher-dimensional spatial geometry for the theoretical description. Systems of that type occur increasingly often, in connection with the use of microelectrodes and their arrays,4 hydrodynamic methods,5 electrochemical sensors,6 or other complicated experimental devices such as, for example, SECM.7 Owing to the computational costs of the CV simulations, there exist a number of important tasks, related to the theoretical modeling and experimental data analysis in CV, that still cannot be performed in a fully satisfactory way. These tasks include the following. (A) Rapid Visualization and Exploration of the Effects of the Model Parameters on the Voltammograms. The voltammograms usually depend on many model parameters, such as reaction rate constants, charge-transfer coefficients, conditional electrode potentials, diffusion coefficients, initial concentrations, etc. Investigating and understanding the effects of these parameters on the transient curves is a necessary element of CV studies and an important component of electrochemical education. Proper visualization could substantially facilitate these tasks, as is the case with treating other data in analytical chemistry.8 However, thus far it has not been possible to create graphical representations of the voltammograms capable of immediately responding to any interactively imposed changes of the parameters. (B) Rapid Estimation of Multiple Parameters from Experimental Voltammograms. Currently employed methods for fitting of the theoretical voltammograms (or other electrochemical responses) to the experimental curves include the following: (3) Britz D. Digital Simulation in Electrochemistry, 3rd ed.; Springer: Berlin, 2005. (4) See, e.g.: Forster R. J. In Encyclopedia of Electrochemistry: Bard A. J., Stratmann M., Unwin P. R., Eds.; Wiley: New York, 2003; Vol. 3, Chapter 2.5, p 160. (5) See, e.g.: Mount A. R. In Encyclopedia of Electrochemistry; Bard A. J., Stratmann M., Unwin P. R., Eds.; Wiley: New York, 2003; Vol. 3, Chapter 2.4, p 134. (6) See, e.g.: Baker E. Anal. Chem. 2004, 76, 3285. (7) See, e.g.: Horrocks B. R. In Encyclopedia of Electrochemistry; Bard A. J., Stratmann M., Unwin P. R., Eds.; Wiley: New York, 2003; Vol. 3, Chapter 3.3, p 444. (8) Wolkenstein M. G.; Hutter H.; Grasserbauer M. Trends Anal. Chem. 1998, 17, 120.
Analytical Chemistry, Vol. 78, No. 6, March 15, 2006 1807
gradient-based methods (e.g., the Gauss-Newton method implemented in the DigiElch9 program), gradient-free methods (e.g., the simplex method implemented in CVFIT10 and ELSIM11 programs), and global search methods (e.g., genetic algorithms implemented in the PIRoDE12 program). All these methods rely on simulating many transient curves, with the aim of gradually approaching a parameter combination that ensures the best agreement between the model solution and experimental data. Thus, time-consuming simulations can result in intractable parameter estimation efforts. (C) Real-Time, On-Line Data Analysis of Experimental Voltammograms. As a consequence of (B), the theoretical analysis of experimental voltammograms is still performed offline (i.e., after the experimental measurements are completed) and over a much longer time scale. For example, model discrimination and identification, which often relies on fitting of a series of conceivable models and choosing the best one, takes even longer time than parameter fitting for a single model. Such practice is not satisfactory, if the development of fully automatic computeraided electrochemical investigations,13 including real-time measurements14 and treatment of the data, is to become a reality. Instead, computer-aided modules for rapid on-line theoretical data analysis should be designed and included into devices such as polarographs, or even (voltammetric) sensors, to enable a variety of relevant chemometric procedures15,16 in the course of the data acquisition. This objective has not been possible to achieve with the current speed of CV simulations. In view of the above needs, it would be desirable to have a nonsimulative method that would allow one to compute CV curves with a high speed, independently of the complexity of kinetic models. Techniques of this kind have already been developed in other areas of chemistry, where computationally expensive simulations occur, for example, in combustion chemistry, atmospheric chemistry, or geochemistry. They are known under the name of solution mapping17-20 and consist in preparing a collection of easily computable algebraic expressions that closely approximate the dependencies of the various model responses (e.g., dynamic concentration profiles and their functionals) on the parameters over a domain of interest. In the present study, we examine the utility of solution mapping for the fast computation of cyclic voltammograms. Various forms of solution mapping have been discussed in the literature,17-20 but the remainder of this paper will focus on solution maps that are look-up tables of certain coefficients corresponding to the nodes of a suitably chosen discrete grid in the parameter (9) http://www.digielch.de. (10) Gosser D. K., Jr. Cyclic Voltammetry, Simulation and Analysis of Reaction Mechanisms; VCH: New York, 1993. (11) Bieniasz, L. K. Comput. Chem. 1997, 21, 1. (12) http://www.elsyca.com. (13) Bieniasz L. K. In Modern Aspects of Electrochemistry, Conway B. E., White R. E., Eds.; Kluwer/Plenum: New York, 2002; Vol. 35, p 135, and references therein. (14) See, e.g.: Brett C. M. A. Electroanalysis 1999, 11, 1013. (15) See, e.g.: Brown S. D.; Bear R. S. Crit. Rev. Anal. Chem. 1993, 24, 99. (16) See, e.g.: Pravdova´ V.; Pravda M.; Guibault G. G. Anal. Lett. 2002, 35, 2389. (17) Dunker A. M. Atmos. Environ. 1986, 20, 479. (18) Frenklach M.; Wang H.; Rabinowitz M. J. Prog. Energy Combust. Sci. 1992, 18, 47. (19) Tura´nyi T. Comput. Chem. 1994, 18, 45. (20) Pope S. B. Combust. Theory Modell. 1997, 1, 41.
1808
Analytical Chemistry, Vol. 78, No. 6, March 15, 2006
space. In the simplest, “brute force” look-up table case, the coefficients represent tabulated values of the simulated responses, with each response corresponding to a different grid node. Approximate responses corresponding to off-grid combinations of the parameter values are then obtained by interpolation. However, other kinds of look-up tables are conceivable, as we shall see below. The construction of any such response map requires a single effort of performing many simulations, but then the map can be stored and reused multiple times without the need to repeat the costly simulations. In the case of time- and space-dependent model responses, an additional tabulation in the time and space domains is likely to be needed, apart from the tabulation in the parameter domain. To avoid complications associated with the interpolation in time and space, it is convenient if the spatiotemporal tabulation grids do not have to vary upon a change of the parameters. Cyclic voltammograms, which are responses dependent on time (or electrode potential), are well suited for such fixed-grid tabulations. They usually do not exhibit any dynamic local narrow structures that would require a time (or potential) tabulation grid modification upon changing the model parameters. Exceptions from this rule exist, but are rare (e.g., current spikes possible with the EE-DISP reaction mechanism21 that are best resolved by adaptive grids22). In contrast, other transient methods such as, for example, potential-step chronoamperometry, often exhibit narrow initial time layers if fast homogeneous reactions are present, so that they may require intensive temporal tabulation grid rearrangements when the homogeneous rate constants vary over large intervals. For similar reasons, finding a fixed grid for the spatial tabulation of the concentration profiles of the chemical species (in any kind of transient experiments) can often be difficult, owing to the frequent presence of narrow boundary and interior layers for fast homogeneous reactions.3 Therefore, although in principle solution mapping could be applied to other transient methods and to spatially dependent concentrations as well, the voltammograms should be the easiest to map. Look-up table solution mapping is not entirely new to electrochemistry: Speiser and co-workers23-26 created tables representing the dependencies of the various features of the voltammograms (such as peak potential, height, and width) on model parameters. Combined with suitable interpolation schemes, the tables were later used for parameter estimation. Also, Alden and Compton27,28 used a similar approach to construct working curves and surfaces (of responses dependent on up to two parameters) for the analysis of steady-state voltammetry at hydrodynamic electrodes. In the present study, we take a further step and consider maps representing entire CV transients, without imposing limitations on the number of model parameters. The main difficulty, inherent in solution mapping, is known as “the curse of dimensionality”, which refers to the possible exponential growth of the volume of the simulated data that needs (21) Dietrich M.; Heinze J.; Krieger C.; Neugebauer F. A. J. Am. Chem. Soc. 1996, 118, 5020. (22) Bieniasz L. K. Electrochem. Commun. 2002, 4, 5. (23) Speiser B. Anal. Chem. 1985, 57, 1390. (24) Scharbert B.; Speiser B. J. Chemom. 1988, 3, 61. (25) Speiser B. J. Electroanal. Chem. 1991, 301, 15. (26) Speiser B. Trends Anal. Chem. 1991, 10, 9. (27) Alden J. A.; Compton R. G. J. Phys. Chem. B 1997, 101, 9741. (28) Alden J. A.; Compton R. G. Electroanalysis 1998, 10, 207.
to be generated and stored, with the number P of model parameters. For example, in the brute force look-up table approach, if one takes N grid nodes along every parameter axis, then NP model responses need to be created and stored. Assuming, for example, that N ) 10 and that every response consists of at least M ) 100 discrete double-precision values represented by four bytes of memory, it is easy to see that, even with a moderate number of about P ) 9 parameters, the size of the table would exceed the typical capacity of the currently available hard disks. The amount of time needed for creating the table would also be unacceptable. Further problems arise with the accuracy of multidimensional interpolation, when the number of parameters increases. These difficulties have appeared insurmountable; however, the recently described “high-dimensional model representation” (HDMR)29-31 offers one practical approach to overcome them. Therefore, in the present work, we will use the HDMRbased look-up tables for solution mapping in CV. The following section explains the principles of the HDMR technique. Further sections contain CV models on which the HDMR solution mapping has been tested, computational details, discussion of the maps obtained, and conclusions. HDMR TECHNIQUE The theoretical foundations of HDMR have been described in considerable detail in several papers29-31 and need not be repeated here. We only summarize the essentials relevant for the present work. Let a (generally M-dimensional vector-valued) function f(p1, p2, ..., pP) represent a particular model response (e.g., a sequence of M tabulated CV current values corresponding to a series of selected time or potential values) as a function of P parameters p1, p2, ..., pP. The HDMR is an expansion of multivariate functions into an additive superposition of functions, where each contains fewer variables. In other words, a multivariate function is expressed as a sum of correlated functions with hierarchically increasing numbers of arguments, up to the maximum number (here P). The term “correlated” used here is to distinguish it from that employed in statistics and means a joint effect of several variables in the multivariate function. Hence, in the present case
f(p1,p2,...,pP) ) f0 +
∑ f (p ) + ∑ i
1eieP
i
fij(pi,pj) +
1ei
