A Comparison of Finite Difference Schemes for Computational Modelling of Biosensors

Abstract. This paper presents a one-dimensional-in-space mathematical model o f an amperometric biosensor. The model is based on the reaction-diffus ion equations containing a non-linear term related to Michaelis-Menten kinetics of the enzym atic reactions. The stated problem is solved numerically by applying the finite diff rence method. Several types of finite difference schemes are used. The nu merical results for the schemes and couple mathematical software packages are compared a nd verified against known analytical solutions. Calculation results are compared in terms of the precision and computation time.


Introduction
The interest in biosensors is constantly growing as the range of practical applications of electrochemistry increases.Biosensors are small analytical devices capable of detecting specific compounds and therefore they are often applied in the fields of clinical, industrial, environmental and agricultural analyses [1][2][3].A biosensor device is composed of biologically responsive material, mostly enzymes, and the electrode.Enzyme interacts with the target substance yielding the product.This process is usually described by Michaelis-Menten kinetics of the enzymatic reactions [4][5][6].Amperometric biosensors are based on the measurement of the Faradaic current when a constant potential is kept.The current on electrode results due to the direct oxidation or reduction of an electroactive spieces.
Analytical solutions for mathematical models of the biosensors are obtainable exclusively in special cases [7,8].In common case the models have to be solved numerically [9,10].Finite difference method is one of the most popular approximation techniques [11].Numerous types of finite difference schemes can be considered for the solution of nonlinear reaction-diffusion systems [4,11].Three major factors must be taken into account when choosing the technique for simulation: the accuracy of the solution, computation time needed to solve the problem and ease of use of the technique.This work is focused on the analysis of several most commonly known finite difference schemes using computer simulation.

Mathematical model
An amperometric biosensor can be considered as an amperometric electrode, having a layer of enzyme immobilized onto the surface of the electrode.We assume the symmetrical geometry of the electrode and homogeneous distribution of the immobilized enzyme in the enzyme membrane.
We consider the following enzyme-catalysed reaction In this scheme the substrate (S) combines reversibly with an enzyme (E) to form a complex (ES).The complex then dissociates into a product (P) and the enzyme is regenerated.
Assuming the quasi steady state approximation, the concentration of the intermediate complex (ES) do not change and may be neglected when simulating the biochemical behaviour of biosensors [1,2].The scheme (1) reduces to a simplified model of enzymecatalyzed reaction, where the enzyme (E) binds to the substrate (S) producing the product (P) is considered, Coupling the enzyme-catalyzed reaction with the one-dimensional-in-space diffusion, described by Fick's second law, leads to the following system of equations [8,9]: where S(x, t) and P (x, t) are the substrate and product concentrations, respectively, t stands for time and x -for space, D S and D P are the diffusion coefficients of the substrate and product, respectively, K M is the Michaelis-Menten constant, V max is the maximal enzymatic rate attainable when the enzyme is fully saturated with substrate, d is the thickness of the enzyme membrane.Let x = 0 represents the electrode surface, while x = d represents the bulk solutionmembrane interface.The biosensor operation starts when some substrate appears on the surface of the enzyme layer, where S 0 stands for the concentration of substrate in the bulk solution.
In the case of amperometric biosensors, due to the electrode polarization, the concentration of the reaction product at the electrode surface is being permanently reduced to zero.The substrate does not react at the electrode surface.If the substrate is wellstirred and in powerful motion, then the diffusion layer (0 < x < d) remains at a constant thickness of d during the biosensor operation.This is used in the boundary conditions given by: The measured current is accepted as a response of an amperometric biosensor in a physical experiment.The current depends upon the flux of the reaction product at the electrode surface, i.e. at the border x = 0. Consequently, the density I(t) of the anodic current at a time t can be obtained explicitly from Faraday's and second Fick's laws using the flux of the product concentration at the surface of the electrode, where n e is a number of electrons, involved in charge transfer at the electrode surface, and F is the Faraday constant.We assume that the system (3)-( 5) approaches a steady state as t → ∞, where I p is the density of the steady state current.

Solution of the problem
The analytical solutions for nonlinear partial differential equations generally do not exist.Equations ( 3)-( 5) describing the action of an amperometric biosensor do not have ones either, so numerical approximation must be used.We applied finite difference technique to solve (3)-( 5) the boundary value problem numerically [10,11].

Analytical solutions
A non-linear term in equations ( 3) turns to linear one in special cases of the substrate concentration, The analytical solutions are known for the boundary value problem ( 3)-( 5) in the cases of linear reaction terms [7,8].Exact solutions are helpful in testing models and assessing accuracy of the solution.
If inequality S 0 K M is satisfied, then the biosensor steady state current can be calculated as follows [8]: where σ 2 is a dimensionless diffusion modulus, Damköhler number, The biosensor response is known to be under mass transport control if the enzymatic reaction in the enzyme layer is faster than the mass transport process [4,8,9].The diffusion modulus essentially compares the rate of enzymatic reaction (V max /K M ) with the diffusion through the enzyme layer (D S /d 2 ).If σ 2  1 then the enzyme kinetics controls the biosensor response.The response is under diffusion control when σ 2  1.At the high concentration of the substrate (S 0 K M ), the biosensor steady state current does not depend on the concentration S 0 of the analyte [7], However, in the intermediate concentration cases, i.e. if S 0 ≈ K M , the analytical solutions are unknown and numerical methods are used to solve the problem [4,9,10].

Finite difference schemes
We introduce an uniform discrete grid ω h × ω τ to simulate the biosensor using finite difference method, where T stands for the duration of the process analysis.The differential equations are discretized in that domain assuming the following definitions: S j i = S(x i , t j ), P j i = P (x i , t j ), I j = I(t j ), i = 0, . . ., N ; j = 0, . . ., M. (14)

Explicit finite difference scheme
Using explicit finite difference scheme for the substrate and product concentrations (3) we obtain the following finite difference equations [9][10][11]: The initial conditions (4) in numerical model has the following form For the boundary conditions ( 5) we obtain: The formulae for calculation of current density (6) becomes thus (0 < j ≤ M ): We consider the density I R of the steady state current calculated at the moment T R We used ε = 10 −5 for the calculations.One of the most important features of the scheme is the stability [11].The prerequisite for the stability of the explicit finite difference scheme (15)-( 17) is the following condition: Because of these stability conditions, a number of the time steps must be magnified strongly as the number of the space steps is increased.This leads to the inefficient calculations.

Implicit 1 finite difference scheme
Mathematical model of a biosensor can be solved using several implicit finite difference schemes [9][10][11].Let us name first of them "Implicit 1 finite difference scheme" and compose its equations.
The model governing equation for the substrate concentration in (3) is replaced by the following finite difference equation: Known concentration values of the substrate at the upper layer can be used for calculation of the product concentration.Hence, governing equation for the product concentration can be approximated with: The rest equations ( 4)-( 6) take the same form as those of the explicit scheme.

Implicit 2 finite difference scheme
The model equation for the substrate concentration in (3) may be approximated with a more implicit scheme than equation (21).At a numerator of reaction term the concentration of substrate can be used in a upper level, Present numerical equation is linear, same as (21).The other equations match those obtained using the explicit scheme.

Hopscotch scheme
Using Hopscotch scheme unknown grid points are obtained at two phases [11].On the first phase, even grid points (S[e] j+1 i , P [e] j+1 i ) are calculated explicitly using known lower layer values (S j i , P j i ), On the second phase, odd grid points ) are calculated using odd values of the lower level and already known even values of the upper layer, Computation continues alternating the calculation order of the odd and even points.Hopscotch scheme is fully explicit yet unconditionally stable for V max = 0 and therefore it can operate with any size of time and space steps.

Mathematical software packages
The major mathematical software packages provide tools for solving systems of the partial differential equations.However, a greater ease of use and wider range of solvable problems often comes at the expense of lower precision or less efficiency.
We used Maple (Maplesoft, Inc.) version 10 general-purpose solver "pdsolve" to find numerical solution for the system of the partial differential equations [13].This solver uses finite difference method and can be configured with eleven classical schemes, calculation step size and other parameters.
The biosensor action was also simulated with MATLAB (The MathWorks, Inc.) software package [14].The problem was solved using built-in solver "pdepe", which provides a numerical solution for systems of differential equations in single spatial dimension and time.

Results and discussion
Computer simulation was used to compare accuracy and performance of the solution techniques.Since the system of linear algebraic equations is tridiagonal it can be solved efficiently [11].Calculation results are compared in terms of precision and computation time.We define the relative error E as the absolute difference of the steady state current density estimated by analytical and numerical solutions divided by the steady state current density of analytical solution.
where I R is the numerical solution defined by (19) while I l and I g are analytical solutions defined by ( 10) and ( 12), respectively.
The following values of the model parameters were employed: The routines of the finite difference method were implemented in Java programming language [12].For performance reasons, we executed programs of the mathematical software packages using command-line approach.The experiments were performed on the 2 GHz Intel Core 2 Duo Processor with 1GB of RAM.
As a first test problem, relative errors of the finite difference schemes and mathematical software packages were examined using two known analytical solutions (10) and (12).We applied M = 10 2 for the calculations using implicit 1, implicit 2 and Crank-Nicolson schemes and M = 10 5 using explicit and Hopscotch schemes because of the stability constraints on the time step.All the considered finite difference schemes yield very similar precision, therefore only explicit scheme results are presented in Fig. 1 as the example.The smallest relative errors are obtained using finite difference schemes.Maple package calculates the steady state biosensor current more accurately than MATLAB.In cases of high substrate concentration (12) Maple's results are as precise as those obtained by explicit scheme.The numbers of steps used in calculations do not influence the accuracy of the MATLAB solution.In Fig. 2 the finite difference schemes are compared to the explicit scheme.The Hopscotch scheme differs very slightly from the explicit scheme in both analytical solution cases, whereas the implicit schemes showed the maximal difference of approximately 0.8 % when the number of space steps N equals to 160 (Fig. 2).The relative errors computed using the analytical solution (12) at S 0 K M are by a few orders of magnitude larger than the corresponding errors at S 0 K M calculated using (10).This could be explained by less accuracy of the analytical solution ( 12) compared to the solution (10) [7,8].Considered schemes yield more similar results using analytical solution (12).
In the next test problem, we consider the computation time as a function of the relative error (Fig. 3).Introducing different limits for the relative error E, the computation time T E ( ) is given by: where T N,M is the time of calculation at given numbers of grid steps N and M .T E ( ) is the minimal time of computation needed to achieve the relative error E not greater than .The calculations were performed for very different values of space and time steps, N, M ∈ {20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 10240}.
As one can see in Fig. 3, the implicit and Hopscotch are the fastest schemes to achieve the required relative error.Despite being very computationally intensive, partial differential equation solvers in mathematical software packages cannot accurately calculate the results.The MATLAB solver does not obtain higher precision than 0.1.
Finally, the computation times for different grid steps are reported in Tables 1, 2. The Table 1 shows that there is a very small difference between schemes implicit 2 and Crank-Nicolson and they are the most computationally intensive schemes.Explicit scheme is the fastest computation technique.Mathematical software packages are significantly more computationally intensive, particularly Maple, see Table 2

Conclusions
In this article, several finite difference schemes were applied for modelling an amperometric biosensor.Using all the considered schemes quite satisfactory results were obtained when sufficient number of steps of the discrete grid is employed.The best accuracy is achieved using implicit calculation and Hopscotch approaches.For the problems where accuracy is not a significant factor but the speed is, the simplest explicit scheme should be used.
General-purpose solvers of Maple and MATLAB are less precise to simulate the biosensor action and need more computation time.Those solvers can be applied for basic problems while taking advantage of the simplicity.

Table 1 .
. Computation time [ms] by the finite difference schemes, N = 100

Table 2 .
Computation time [s] by the mathematical software packages, N = 100