Relaxation and Simplex mathematical algorithms applied to the study of steady-state electrochemical responses of immobilized enzyme biosensors: Comparison with experiments

A description of the implementation of the relaxation method with automatic mesh point allocation for immobilized enzyme electrodes is presented. The advantages of this method for the solution of coupled reaction–diffusion problems are discussed. The relaxation numerical simulation technique is combined with the Simplex fitting algorithm to extract kinetic parameters from experimental data. The results of the simulations are compared to experimental data from self-assembled multilayered electrodes comprised of glucose oxidase (GOx) and an Os modified redox mediator and found to be in excellent agreement.


Introduction
For amperometric enzyme electrodes with the enzyme entrapped in a film, the interplay between diffusion and kinetics results in highly non-linear differential equations for which there are no closed form analytical solutions. Approximate analytical solutions for selected limiting cases have been derived [1,2]. However, it is very useful to complement these with fast and reliable numerical simulations which can treat the intermediate cases.
The main advantages of numerical simulation are that it can be used to calculate the amperometric response over the whole range of experimental parameters and that it gives calculated concentration profiles for the different species within the film. In these respects it complements the approximate analytical treatments available in the literature that are only valid under certain limiting conditions, i.e. for limited ranges of the experimental variables. Numerical techniques are particularly valuable for exploring the behaviour of the system in the regions between the analytical limiting cases where, typically, two parameters for the system are comparable and a limiting behaviour is not well established.
Ideally the most powerful approach is to combine the use of approximate analytical solutions that provide physical insight into the nature of the rate limiting processes and the physical behaviour of the system with numerical approaches that can be used to fit the full range of experimental data and to extract the best estimates of the controlling kinetic parameters. For amperometric enzyme electrodes this offers the possibility of rational electrode design where the model is used to predict the amperometric response of the electrode [3]. It also offers the possibility to investigate detailed questions about the operation of the electrode such as the true activity of the immobilised enzyme or the direct determination of the rate of mediated electron transfer to the enzyme within the film -issues which are highly relevant to the optimal design of biosensors and efficient biofuel cell electrodes. Bartlett and Pratt [1] used numerical simulations to validate their approximate analytical solutions and to investigate the boundary regions between the different limiting cases in a clear example of how numerical simulations and approximate analytical techniques can complement each other in the modelling of a particular physico-chemical system.
Fitting experimental data to theory, both approximate analytical treatments and numerical simulations, is essential both to test the validity of the model and to extract the relevant kinetic constants such as the various rate, diffusion and partition coefficients. This is most effectively achieved using non-linear least squares optimization to fit the experimental data to the theory. To achieve this in a reasonable time using a numerical model requires a computationally efficient numerical simulation algorithm since the simulations will have to be repeated many times during the iterative fitting of the kinetic parameters.
Several numerical techniques have been used to solve both the transient and the steady-state situations with various boundary conditions. For a survey of the early literature see [3]. More recently, the problem of amperometric enzyme electrodes, has been analysed by digital simulation in several papers [1,[4][5][6][7][8][9][10][11][12][13][14][15]. As far as we are aware, Bartlett and Pratt [1] are the only ones to have used the relaxation method [16] to treat an electrochemical problem. The advantage of the relaxation method is that it gives a rapid and stable solution for the steady-state behaviour of strongly coupled non-linear differential equations. In their earlier work Bartlett and Pratt did not describe the implementation of the technique in detail. However, they found excellent agreement between their simulations and the calculations from their approximate analytical treatment although they did not compare their simulations with experimental.
In this paper, we build on and significantly extend the work of Bartlett and Pratt. We describe in detail the relaxation method and its implementation to solve the coupled non-linear kinetic-diffusion problem for an amperometric enzyme electrode. We show how the relaxation method can be combined with the Simplex method [16], in order to fit experimental data and extract the relevant kinetic parameters. We show that the simulations are in good agreement with experimental data from electrostatically self-assembled multilayers of Glucose Oxidase (GOx) and an Os modified polyelectrolyte (PAH-Os), and we use our combined relaxation simulation and Simplex optimisation to determine kinetic data for the system. This new treatment should be applicable to enzymes immobilized in modified polymers, in redox polymer films [17], in redox hydrogels [18], by antigen-antibody interaction (for example avidin-biotin) [19][20][21], and electrostatically self-assembled systems in general [22,23]. The treatment is equally valid for both thin monolayers as well as thick multilayered films and offers an efficient approach to the treatment of experimental data. The analysis presented here, together with appropriate approximate analytical formulae can be used in the design and optimization of enzyme electrodes for analytical and biofuel cell applications when the final objective is to obtain a detailed understanding of the operation of the electrode and to be able to predict the response.

The model
A general kinetic model for an enzyme-membrane electrode has been described previously [1] and is briefly reviewed in Fig. 1. For fuller details the reader is referred to the original paper. The model describes the general situation of a redox mediated immobilized enzyme electrode. In principle the redox mediator could be co-immobilized in the film or it could be present in the bulk solution and undergo partition into the film. In the latter case, the mediator could be present in the bulk solution either in its oxi- and substrate (S) occurs within the film with diffusion coefficients D e À (electron hopping) and D s respectively. Partition of substrate between the film and the bulk solution is described by the partition coefficient K s . The homogeneous enzyme kinetics occurs throughout the film from v = 0 to v = 1. The reduced mediator is reoxidized to produce A at the electrode surface. E 1 and E 2 are the oxidized and reduced enzyme respectively. dized or reduced form. In the present paper we only consider the situation where the mediator is immobilized within the film. This is a very common situation exemplified by, for example, the popular and successful approach of using a redox hydrogel pioneered by Heller and his colleagues [24]. The extension of our treatment to either of the other two situations can be achieved by an analogous numerical treatment taking account of the appropriate boundary conditions [3].
The reactions occurring in the film can be written and at the electrode where A and B are the oxidized and reduced forms of mediator, E 1 and E 2 are the oxidized and reduced forms of enzyme, S and P are the substrate and product of the enzymatic reaction and ES is the complex formed between enzyme and substrate (see Appendix 1 for a full list of the symbols used in this work). The substrate undergoes partition between the bulk solution and the film (partition coefficient K S ) and then diffuses within the film with a diffusion coefficient D S . The mediator is assumed to be confined within the film. Either it can physically diffuse within the film described by a diffusion coefficient D A or, if the mediator is covalently bound within the film, charge propagation occurs by electron hopping self-exchange between the reduced and oxidized dorms of the mediators described by a diffusion coefficient D e . According to the Dahms-Ruff formalism [25,26] these two situations are equivalent with D A = D e . Michaelis Menten kinetics are assumed for the enzyme-substrate reaction. The oxidized mediator regenerates the enzyme with a mediator-enzyme reoxidation constant k, according to the conventional ''ping-pong" mechanism [27].
The second-order differential equations describing the system in the steady-state are The symbols in brackets refer to the concentrations of the corresponding species. These concentrations vary the position within film. Eqs. (4) and (5) are non-linear second order differential equations and have no closed form analytical solution. In this paper we solve these equations using the relaxation method. First we recast the problem in terms of the following dimensionless variables [1,3]

Numerical simulations
To analyse the steady-state enzyme electrode response we need only simulate the steady-state solution to Eqs. (4) and (5). The high rates of reaction, and thus steep concentration gradients which can occur within the immobilized layer at the electrode surface, make simulations using an explicit finite difference method impractical. Implicit or semi-implicit methods suffer the drawback that the complex kinetic scheme cannot be solved directly. Consequently we chose to use a steady-state simulation method. The obvious choice, the shooting method [16], is inappropriate in the present case because of the combination of boundary conditions at the two membrane interfaces -the electrode surface and the outer surface of the film. On the other hand, the relaxation method [16] can simultaneously solve any number of coupled first order differential equations, provided that sufficient boundary conditions are known. Eqs. (4) and (5) contain second order diffusional terms but this is not a problem since second (or higher) order differential equations can be expressed as a combination of first order terms.
In the relaxation method a set of ordinary differential equations is replaced by a set of approximate finite difference equations on a grid of points which spans the domain of interest, in this case, the enzyme layer. The principle behind the relaxation method is to start with an initial guess for the concentration profiles within the enzyme layer. The method then computes the resulting errors at each grid point and uses these computed errors to make an improved guess. Hence the name since it could be said that the algorithm relaxes to the correct solution. However it is important to note that these intermediate solutions, generated en route to the final steady-state solution have no physical significance and that the final solution is independent of the initial guess.
To illustrate the method consider the trivial case of a single equation W in terms of a variable y We starting with an initial guess for the value of y and this is adjusted iteratively by a small value Dy. For the correct solution Alternatively this can be written Rearranging Eq. (16) therefore gives the value of Dy in terms of the previous estimate for y: When dealing with a large number of coupled equations and variables V is a matrix containing the differentials of all equations W with respect to each variable y. Solution of the equations therefore requires inversion of this matrix. Further details are given by Press et al. [16]. When performing electrochemical simulations using finite difference methods a useful reduction in the number of points required for a given accuracy (and hence increase in the speed of the simulation) can be achieved by using a non-linear mesh spacing [28]. This concentrates the computational effort in those regions where the concentrations and fluxes change most rapidly. This is normally achieved by using a high mesh density, i.e. closely spaced points, where the concentration gradient is steepest. For simple electrochemical systems such as the ec or the ec 0 mechanisms, or where a rotating disc electrode is used, the concentration profile is steepest at the electrode surface. Thus Feldberg [29] used a mesh spacing based on an error function to give a close point spacing at the electrode, with the spacing increasing with distance from the electrode in an appropriate way. This approach, and similar methods, are widely used. The problem with this type of pre-defined mesh spacing is that it requires an a priori knowledge of the approximate shape of the concentration profiles. In the present case [1], the concentration profiles may be steep at the electrode surface, at the membrane-solution interface, at both, or even at some point in between depending on the balance between the different rate processes. It has also been found that a high mesh density is required where the concentration gradients are rapidly changing, i.e. where d 2 [A]/dx 2 or d 2 [S]/dx 2 are large. For concentration profiles such as those shown if Fig. 2, this requires a higher mesh density at some region within the film that is not necessary near the boundaries. The relaxation method allows the optimum mesh spacing to be achieved automatically [16] by deriving an equation describing the relationship between the mesh density and the concentration profiles. This equation is then solved simultaneously along with Eqs. (4) and (5). Fig. 3, as described in the legend, a linear mesh spacing would not accurately describe the profile, since the reaction layer only covers approximately 1% of the film, near the electrode surface. A pre-defined mesh spacing with, for example, mesh density decreasing exponentially with v, would be better than nothing, but, since the substrate profile may be steep near v = 1, a high mesh density would also be required there. This would be wasteful of points, and computationally inefficient, if the concentration profiles are only steep at one of the boundaries. Fig. 4 shows the effect of making the mesh density depend on the second differential of the concentration profiles, i.e. the rate of change of the gradient. Here a high mesh density is required at a position away from the boundaries. Since this may occur anywhere within the film, a single pre-defined mesh spacing function would not be useful.

Figs. 3 and 4 demonstrate the benefits of automated mesh point allocation. In
Both the relaxation and Simplex methods and their implementation are described in detail by Press et al. [16]. In this paper we only describe how these methods are used to simulate the immobilized enzyme electrode and to analyse experimental data. This should be sufficient to allow the reader to adapt the program for other model systems.
Our programs make use of the subroutines published by Press et al. [16] and are available from the authors on request.
3.1. Implementation of the relaxation method 3.1.1. Representation of concentration profiles In our system, four variables are required to describe the concentration profiles of mediator and substrate here only a small number of points, i.e. a low mesh density, is sufficient. In region B, the gradient is high, so a higher mesh density is required. In region C, the gradient of the concentration profile is changing rapidly. This also requires a high mesh density.
Then d 2 a dv 2 is represented as the first order differential equation dy 3 dv .

Automated mesh point allocation
The relaxation method uses a grid of N mesh points to represent the distance 0 6 v 6 1. Each variable y r above is therefore an array of values where 0 6 n 6 N.
In order to implement an adaptive non-linear mesh spacing we require three more y r,n arrays: where q is equal to the point number n, i.e. it ranges from 1 to N. The variable Q is proportional to q, but does not have a defined range of values. The relationship between Q and distance v at any one point is determined by the . Graph (i) shows the mesh density, dq/dv (=h/Dy 7 ). The mesh density is higher near to the inner boundary than it is near to the outer boundary, due to the relative gradients of a and s. The peak in the mesh density at v = 0.017 is due to P 4 and P 5 in Eq. (25). This gives a high mesh density where the gradients of the concentration profiles change rapidly, i.e. at the titration point. The effect of this is shown by graph (iii), which shows the region of the concentration profile around the titration point. The distance 0.16 6 v 6 0.18 contains 8 points, compared with an average of approximately 4 points over the same distance in the rest of the film. Simulated for j = 200, c = 0.2, g = 1, l = 0.001, a e = 1, with weighting parameters, P 1 = 50, P 2 = 2.5, P 3 = 2.5, P 4 = 100, P 5 = 100. mesh density function U ¼ dQ dv . A large value of U indicates that a high mesh density is required, i.e. a large number of points within a given distance Dv. The definition of U given by Press et al. [16] is not suitable for concentration profiles, so a new definition of U was used. A high mesh density is required where the concentration profiles are steep, U should therefore be proportional to da/dv and ds/dv, so that it is large in, for example, region B of Fig. 2.
For our system, it was also found to be necessary to have a high mesh density where the gradient of the profile was rapidly changing, i.e. in region C of Fig. 2. If we analyse what Bartlett and Pratt [1] call the titration case, i.e. a situation where both mediator and substrate are consumed within the film resulting in a change from mediator to substrate limited kinetics through the film, there is a very sharp change in gradient at the titration point, v * .
Taking account of these considerations, we find that a suitable definition of U is which, in terms of the program variables y, becomes: The variables P 1 to P 5 are weighting parameters, the larger the value the more effect that part of the equation has on the mesh spacing. The P variables do not need to be normalised, since the value of Q is automatically adjusted so that as 1 6 q 6 N then 0 6 v 6 1. Since it is the magnitude, and not the sign, of the differentials in Eq. (25) which determines the mesh spacing, absolute values are used. With reference to Fig. 2, a high (relative) value of P 1 favours a linear point spacing. A high value for P 2 favours a high mesh density in region B. A high value of P 4 favours a high mesh density in region C. P 3 and P 5 are the equivalents to P 2 and P 4 but for the substrate concentration profile.
The reason for using the method shown in Eqs. (25) and (26) instead of d 2 a/dv 2 and d 2 s/dv 2 for the second differentials, is that the latter would give the following expression in which the expression for the mesh density, 1/Dy 7 (=1/ Dv) occurs in the definition of the mesh density function, U. Errors in y 7 would thus cause divergent behaviour.
The form of U shown in Eqs. (25) and (26) should be generally suitable for any system of concentration profiles.

Differential equations
As described by Press et al. [16], G variables require G differential equations W to solve them. Two of these are given by the relationships between y 1 and y 3 and y 2 and y 4 similarly; Two more are given by Eqs. (4) and (5): From Eq. (22) As described earlier, there is a linear relationship between q and Q. Therefore The final differential equation is given by the mesh spacing function U 3.1.4. Finite difference representation These seven equations are expressed in finite difference form as W g,n , which couples two points n and n À 1, such that dy r;nÀ1=2 dq % y r;n À y r;nÀ1 ¼ y m r ð35Þ and y r;nÀ1=2 % y r;n þ y r;nÀ1 The abbreviations y m r and y p r are used to simplify the notation.
Using the notation shown to the right of Eqs. (35) and (36), the seven finite difference equations are therefore: W 7;n ¼ y m 7 À y p 6 h 2ðP 1 þ 1=2P 2 jy p 3 j þ 1=2P 3 jy p 4 j þ P 4 jy m 3 j þ P 5 jy m 4 jÞ ð43Þ where each W g,n should be equal to zero for the correct solution. The variable h is equal to 1/(N À 1), i.e. it is the average point spacing Dv.

Solution of the equations
For the correct solution, all W g,n should be zero. Taking all W g,n , y r,n , and y r,nÀ1 into consideration, Eqs. (15) and (16) become [16]: W g;n ðy r;n þ Dy r;n ; y r;nÀ1 þ Dy r;nÀ1 Þ % W g;n ðy r;n ; y r;nÀ1 Þ þ X G r¼1 ðV g;r Dy r;nÀ1 þ V g;rþG Dy r;n Þ ð44Þ where V g;r ¼ oW g;n oy r;nÀ1 ð45Þ V g;rþG ¼ oW g;n oy r;n ð46Þ in which 1 6 r 6 G, 1 6 g 6 G, 1 6 n 6 N. Evaluation of all Dy r,n and Dy r,nÀ1 therefore requires inversion of the V matrix. This is performed using subroutines published by Press et al. [16] and will not be discussed further here.
It can be seen from Eqs. (45) and (46) that the simulation requires prior evaluation of the differentials of every V g,n with respect to every y r,n and y r,nÀ1 . Thus with seven equations and seven variables, 98 differential equations are required. A few of these will be shown here as examples. For equation W 1 : For r = 1 we differentiate Eq. (37) with respect to variable y 1 at point n À 1 W 1;n ¼ y m 1 À ðy m 7 y p 3 Þ=2 ð48Þ V 1;1 ¼ dðy 1;n À y 1;nÀ1 À ðy m 7 y p 3 Þ=2Þ dy 1;nÀ1 ð49Þ For r = 2 For r = 3 V 1;3 ¼ dðy m 1 À ðy m 7 ðy 3;n þ y 3;nÀ1 ÞÞ=2Þ dy 3;nÀ1 ð53Þ The remaining four V 1,r differentials are evaluated similarly. V 1,r+1 to V 1,r+7 are differentiated with respect to y 1 at point n V g;rþG ¼ oW g;j oy r;j ; 8 6 r þ G 6 14 ð55Þ For example, V 1,14 is W 1,n differentiated with respect to y 7,n V 1;14 ¼ dðy m 1 À ððy 7;n þ y 7;nÀ1 Þy p 3 Þ=2Þ dy 7;n ð56Þ Thus we can see that many of the 98 differential equations are either 0, or yield a fairly simple expression so that the task is not so arduous as it might first appear. A few of the equations are slightly more complicated, for example V 3,1 is Application of the quotient rule gives 3.1.6. Differentiation of absolute values When determining V 7,3 , V 7,4 , V 7,10 and V 7,11 it is necessary to differentiate absolute values (recall the expression for W 7,n Eq. (43)). The strict mathematical way to do this is as follows where sign(f(x)) is +1 for f(x) > 0, À1 for f(x) < 0 and 0 if f(x) = 0. However, in a limited number of cases this caused numerical problems when it was implemented. After extensive testing we found that by taking we were able to obtain robust results for all cases (note: the use of Eq. (61) only affects the automated grid point allocation and it does not introduce any errors into the simulated concentration profiles or current, it only alters the efficiency of the simulation). Results obtained in this way are in good agreement with the analytical solutions corresponding to the different limiting cases. The equation used for V 7,3 is therefore This deals with the points between the electrode and the membrane-solution interface, 2 6 n 6 N À 1. We next consider the boundaries.

Boundary conditions
At the boundaries, the technique is similar. Since there is no n = 0 point, a differential equation for y r,nÀ1 is not required at the inner boundary where n = 1. Similarly, as there is no N + 1 point, the equation W g,n is only differentiated with respect to y r,N . Therefore at the inner boundary v = 0, n = 1: while at the outer boundary v = 1, n = N: where 1 6 r 6 G. Thus there are potentially a total of 2G boundary conditions. To solve G equations, G boundary conditions must be known. There does not necessarily need to be one boundary condition specific to each y array, i.e. some y arrays may have both boundary conditions known, while others have none. For our problem, four boundary conditions are known for y 1 to y 4 . At the inner boundary v = 0, n = 1 a ¼ a e ; W 7;1 ¼ y 1;1 À a e ; V 7;8 ¼ 1 ð65Þ At the outer boundary v = 1, n = N: The substrate concentration has its bulk value Three more boundary conditions are required, these relate to the mesh spacing. At the inner boundary At the outer boundary All other V g,r are zero. Note that the ordering of the index g in the equations W g,n is arbitrary. The ordering of the V values within the matrix is however vitally important. Due to the method used to invert the matrix, the first 4 Â 4 block is used for 'pivoting' and must not be singular for the reasons described by Press et al. [16]. What this means is that each of the first four columns (1 6 r À G 6 4) of the matrix V g,r must contain at least one non-zero value within its four rows (4 6 g 6 7). This is shown diagrammatically in Fig. 5 and is explained in detail by Press et al. [16].

Fluxes
Once the correct values of the y arrays have been found by the simulation, the values of the dimensionless mediator and substrate fluxes are given by where J A is the dimensionless form of the flux measured at the electrode.

Simplex algorithm
The Simplex algorithm [16,30] was used to determine the parameter values which gave the best fit of the theoretical function to the experimental data. The algorithm calculates the goodness of fit using the equation of minimum squares where X is the number of experimental data points and then seeks to minimize f 2 . This fitting algorithm requires a set of starting values for the parameters to be fitted. We generate these by fitting of our experimental data to the approximate analytical expressions [1] that describe the system. In a subsequent paper, we will present a full experimental data set and give a detailed description of the concerted approach used to analyse all of the data combining the numerical simulation methods described here with fitting to the approximate analytical expressions. . All other elements are equal to zero. For pivoting to work, each of the first four columns (1 6 r À G 6 4) must contain at least one non-zero element within its first four rows (4 6 g 6 7) [16]. The matrix (i) shows the initial arrangement of the columns. It can be seen that columns 2 and 3 contain only zeros in their first four rows. Swapping column 2 with 5, and column 3 with 7 corrects this, as shown in matrix (ii).

Incorrect solutions
One potential problem with the relaxation method is that there may be more than one mathematical solution which satisfies the set of equations and boundary conditions. With a transient simulation method, since it proceeds stepwise from the initial conditions, it mimics the physical behaviour of the system. Such a method, correctly implemented, should therefore reach the correct physical solution to the system. Steady state simulations, however, do not follow the physical behaviour of the system, and will converge on the first mathematically correct, or approximately correct, solution which they encounter. Usually if mathematically correct solutions other than the required physical solution exist, these will not correspond to physically possible solution, for example they may give negative concentrations or normalized concentrations greater than 1. Such solutions can therefore be rejected, and the simulation repeated using different starting conditions in order to look for the correct solution. Routines were implemented in the simulation program to avoid such impossible solutions from occurring. A subroutine operates on y 1,n and y 2,n and sets them equal to zero if they become negative, or sets them to 1 if they exceed 1. A second subroutine ensures that the distance must always increase from the electrode surface to the membrane-solution interface. This prevents Dy 7,n from becoming negative during the automated mesh point allocation, this is achieved by checking that y 7,n is greater than y 7,nÀ1 , if this is not true then y 7,n is set equal to y 7,nÀ1 . This results in a step in the y 7 array, from which the simulation can recover to give a smooth y 7 .
For the Simplex method, there also may be more than one set of parameters that minimizes f 2 , i.e. there can be various local minima. It is therefore important to apply sensible constraints on the values of the fitting parameters. These constrains usually derive from three different origins: (1) experimental knowledge about a particular system; (2) approximate values obtained using approximate analytical solutions to which the experimental results can also be fitted; (3) comparison of the fitting parameters obtained when different experimental variables are independently modified for a given system. Therefore, the model is not only simulating the experimental response but also predicting its evolution.
In practice, it is sometimes very difficult to find physically meaningful results when trying to fit three or more parameters at the same time. In these cases it is usually more convenient to start by fixing, or tightly constraining, one parameter until approximate fitting results have been found for the more uncertain ones.
Finally, it is essential to critically evaluate the progress and results of the fitting to check that the resulting parameters are consistent and physically sensible.

Accuracy and limitations of the simulations
The default values for the weighting parameters used in most of our simulations (particularly in all the simulations shown in this paper) are P 1 = 50, P 2 = P 3 = 2.5 and P 4 = P 5 = 100. With the correction subroutines described above and automated mesh point allocation, the simulation works reliably up to j % 1000, for c/g < 1, or j/g 1/2 % 1000 for c/g > 1. This combination of parameters covers all of the physically reasonable situations that could occur for our electrodes. The higher values of j correspond to films that are much thicker than any we are able to achieve experimentally. Over this range of j, g, and c the simulation was in excellent agreement with the approximate analytical solutions with the largest deviations found at the case boundaries where the approximate analytical solutions are least accurate.
Unlike transient simulation methods, the accuracy of the relaxation method does not suffer when high rates of reaction, i.e. large j or j/g 1/2 are used. The simulation simply reaches a point at which it can no longer converge.
Calculations were carried out on a Celeron Ò CPU 2.80 GHz, 448 MB of RAM personal computer using a program written in FORTRAN 77 and incorporating the algorithms for the relaxation and Simplex methods given by Press et al. [16]. For our problem, with a set of seven ordinary differential equations and 200 grid points, each iteration takes less than 1 s and typically less than 100 iterations are required to achieve an acceptable solution.
Gold coated silicon (1 0 0) substrates were employed as electrodes with a 20 nm titanium and a 20 nm palladium adhesion layer and a 200 nm gold layer thermally evaporated using an Edwards Auto 306 vacuum coating system, at p < 1 Â 10 À5 mbar. The freshly evaporated gold film substrates were used once only. To check the quality of the gold surface, the electrodes were cycled in 2 M sulphuric acid between 0 and 1.6 V at 0.1 V s À1 .
Surface modification. An automatic dipping method (Microm DS50 programmable slide stainner from Zeiss Inc.) was used to implement the process described by Hodak et al. [23] to build up layer-by-layer supramolecular multilayers comprised of GOx and PAH-Os. First, the gold surface was modified with sulphonate groups by immersion in a freshly prepared 20 mM 3-mercapto-1-propane sulphonic acid (MPS) solution for 30 min followed by rinsing with Milli-Q Ò water. After thiol adsorption, the first PAH-Os layer was formed by immersion of the thiol-modified Au substrate in a PAH-Os solution for 10 min. The next and subsequent layers were deposited onto the modified surface by alternated immersions in a 1 lM GOx solution and PAH-Os solution respectively for 10 min each, thoroughly rinsing in Milli-Q Ò water at the end of each adsorption step.
A standard three-electrode electrochemical cell was employed with an operational amplifier potentiostat (TEQ-Argentina). A Ag/AgCl; 3 M KCl reference electrode was used together with a large area platinum gauze counter electrode. All electrochemical experiments were carried out in deoxygenated 0.1 M TRIS buffer solutions (0.2 M ionic strength) of pH 7.4, at room temperature (25 ± 2)°C. A SENTECH (Berlin, Germany) variable angle rotating-analyzer automatic ellipsometer (vertical type, 2000 FT model) equipped with a 632 nm laser as polarized light source was employed to measure the film thickness of the electrodes.

Comparison of the simulations to experimental data
We chose to test our numerical simulations on electrostatically self assembled enzyme electrodes. The GOx/ PAH-Os system [22] is a previously well characterized one. It has been extensively studied by cyclic voltammetry, quartz-crystal microbalance, ellipsometry, FT-IR and Raman spectroscopy [32][33][34][35][36]. We know that film thickness, Os surface concentration and enzyme loading all increase with the number of adsorption steps and the catalytic current varies with the film thickness. It has been established that the redox charge propagation within the film is by electron hopping and the diffusion coefficient has been estimated [37]. We know that we can approximate the substrate partition coefficient between the solution and the film, K s , to unity. Finally, we know that, because of the high water content of the films, the glucose diffusion coefficient within the film is almost the same as in pure water [35]. The great advantage of electrostatically self-assembled systems as compared to other enzyme electrodes, such as hydrogels, is that we can vary the design of our electrodes at will choosing from a more or less wide spectrum the thickness, enzyme loading and Os concentration; and, to a lesser extent, k the enzyme-mediator reoxidation rate constant. In this way, we are able to test the simulations on electrodes covering a wide range of parameters, and not only on one or two specific cases.
There are five adjustable parameters, j, c, g, l, and a e , in our model (the values P 1 -P 5 which control the automated grid point allocation are held fixed). Of these five, typically only two or three will be well determined by the curve fitting for any given set of experimental data from a single experiment. Which two or three parameters this is will depend on the case in which the experiment falls. To determine all five parameters with acceptable accuracy the results from more than one experiment carried out over a range of conditions such as film thickness, enzyme loading, mediator concentra-tion corresponding to different cases would need to be combined. We will return to this point in a subsequent paper. Fig. 6 compares the results of the numerical simulations for two different sets of experimental data for the amperometric response as a function of substrate concentration for two self-assembled electrodes of different thickness. Table 1 summarises the numerical results.
Curve a is data for a MPS/(PAH-Os/GOx) 3 electrode. The film thickness measured by ellipsometry was 141 nm. By analyzing the concentration profiles (see Supporting Information) and the ratio of the kinetic constants to the diffusion coefficients, we know that the data points can be fitted to the approximate analytical solution for the thin layer model [1,23] i ¼ where i is the current density, F is the Faraday, C we is the wired enzyme surface concentration and the other parameters have been defined above. Fitting to the thin layer model has been used before to extract k and C we from calibration plots.
This set of data is particularly useful to test the numerical simulations since we can compare the parameters resulting from the combined relaxation/Simplex fitting to the parameters obtained from fitting to the approximate analytical formula that has been previously validated experimentally [23].
When appropriately constrained, the relaxation/Simplex fitting gives the same fitting parameters as a non-linear fitting to Eq. (75). We can see from Fig. 6, that the simulations match our experimental data.
Having tested the approach on a data set for which there is a good analytical approximation we now go on to use it to analyse a data set for which the approximate analytical approach does not work. Curve 'b' in Fig. 6 depicts another set of experimental data for a MPS/(PAH-Os/ GOx) 4 electrode. In this case, the film thickness measured by ellipsometry was 222 nm. When performing the same analysis as for the data in curve 'a' (see Supporting Information), we realize that we cannot fit this data to the thin layer model or any other approximate analytical equation. This is one particular case where the numerical simulations could be of great help in trying to find the unknown parameters for our experimental system. Again from Fig. 6 we can see that there is good agreement between our experimental data and our simulations. If we compare the parameters obtained for the two electrodes (shown in Table 1), we see that [E wired ] for the two films differs by about 20%, but that the two k values are quite different. This can be explained if we consider that in electrostatically assembled multilayer films, the first layers are more strongly affected by the substrate so that thinner films may be expected to show different behaviour from their thicker analogs.

Conclusions
We have proposed, and described in detail, a new numerical simulation method to simulating the steady-state coupled reaction-diffusion problem for immobilized enzyme electrodes. We have investigated the advantages of the relaxation method with non-linear automated mesh spacing over more commonly applied numerical techniques for problems of this type.
The relaxation method is a very powerful tool that could be applied to a wide variety of electrochemical systems, provided that some form of bounded system can be set up, either in the form of a membrane or by enhanced mass transport such as convection or use of a microelectrode.
We have combined the relaxation simulation technique with the Simplex fitting algorithm and used the fitting routine to simulate the amperometric response of self-assembled GOx/PAH-Os electrodes to extract the unknown kinetic parameters from experimental data for the steady-state current at different glucose concentrations. The simulated currents are shown to be in excellent agreement with the experimental results. One of the data sets allowed us to compare the fitted parameters from our relaxation/Simplex fitting routine with parameters obtained by fitting the same data to an approximate analytical formula for the thin layer model. The two methods where shown to be in good agreement.
For the second set of experimental data, there was no unique analytical formula valid for all our data points. In this case, the fitting routine proved useful for finding unknown parameters that characterize our system.