Analysis of Mathematical Modelling on Potentiometric Biosensors

A mathematical model of potentiometric enzyme electrodes for a nonsteady condition has been developed. The model is based on the system of two coupled nonlinear time-dependent reaction diffusion equations for Michaelis-Menten formalism that describes the concentrations of substrate and product within the enzymatic layer. Analytical expressions for the concentration of substrate and product and the corresponding flux response have been derived for all values of parameters using the new homotopy perturbation method. Furthermore, the complex inversion formula is employed in this work to solve the boundary value problem. The analytical solutions obtained allow a full description of the response curves for only two kinetic parameters (unsaturation/saturation parameter and reaction/diffusion parameter). Theoretical descriptions are given for the two limiting cases (zero and first order kinetics) and relatively simple approaches for general cases are presented. All the analytical results are compared with simulation results using Scilab/Matlab program. The numerical results agree with the appropriate theories.


Introduction
Electrochemical biosensors are used as detectors in several commercial analyzers for the accurate and rapid determination of various metabolites such as urea, glucose, lactate, and creatinine [1][2][3][4][5]. These biosensors are fabricated by immobilizing appropriate bioreagents (i.e., enzymes) in a layer adjacent to the sensing surface of the basic electrochemical transducers. The enzyme layer catalyzes the conversion of metabolite molecules, ultimately consuming or producing an electrochemically detectable species. Thus, the analytical performance of these biosensor systems is largely dependent upon the properties of the immobilized enzyme layer incorporated.
A biosensor is an analytical device comprising of a biological element capable to recognize an analyte, coupled with a transducer which generates a signal proportional to the concentration of the analyte and it combines the selectivity and specificity of an immobilized biologically active compound with a signal transducer [6][7][8]. The analytical application of biosensors has become a focus of interest and a subject of rapid progress [9][10][11]. The theory of enzymebased potentiometric sensors is being treated in a series of pioneering contributions. Explicit solutions were derived by Blaedel et al. [12] for the steady state response of such electrodes which apply either to very high or to sufficiently low substrate concentrations. Carr followed the same results from a Fourier analysis as limiting cases for long periods [13,14].
A relatively simple approach was presented by Morf and he obtained an explicit result for the electrode response that applies to the whole range of substrate concentration [15,16]. The response behaviors of potentiometric enzyme electrodes as well as the product release from enzyme reactors were treated for the steady-state case in paper [17] and for the nonsteady state in paper [18] by the same author. Numerical simulations were developed for modeling the reaction and diffusion processes that arise in the functional enzyme membranes of such systems. These simulations represent a kind of virtual experiments and they allowed it to get insight into 2 ISRN Biochemistry the concentration profiles and fluxes of substrate and product species and to analyze the final response characteristics of enzyme-based sensors and reactors [11,19].
To our knowledge, no general analytical expressions of the concentrations of the substrate, product, and current have been reported for all values of parameters [13]. The purpose of this communication is to derive the concentrations of the substrate, product, and current for all values of reaction parameters using homotopy perturbation method. The theoretical treatments make use of the homotopy perturbation method and lead to relationships for all decisive quantities as a function of time. The theoretical results agree with simulated data and offer the basis for reliable predictions of response time ranges for enzyme electrodes and enzyme reactors.

Mathematical Formation of the Problem
In this paper, we consider the analytical system based on an enzyme-containing bulk membrane of thickness that contains a uniform total concentration of the enzyme which is contacted on one side with an aqueous solution of the substrate . The substrate molecules diffuse into the membrane phase where they react in accordance with Michaelis-Menten type enzyme catalyzed reaction [20,21] to yield an electroactive product . Consider where where is the intermediate enzyme-substrate complex, is the number of product species obtained per substrate molecule, 1 , 2 , and 3 are the rate constants of the respective partial reactions, and is the Michaelis constant defined in (2). The influences of reaction and diffusion processes for the species and in the enzyme membrane are described by the following nonlinear governing equations: where [ ] em and [ ] em are the concentrations of the species in the enzyme membrane, and are the corresponding diffusion coefficients, [ ] tot is the total concentration of free enzymes and enzyme-substrate complexes that is assumed to be constant within the membrane including surface zones, is the number of product species obtained per substrate molecule, and 3 is the rate constant for the irreversible step of product formation [17,18]. Now, (3) are solved by assuming the zero fluxes at = 0 and of equilibrium distribution at = [17,18]. The initial state is given by zero concentrations of substrate and product species throughout [17]. Consider For enzyme reactors, the outward flux of the product species at = is described by where is the diffusion coefficient of the product.

Dimensionless Form of the Problem
To compare the analytical results with the simulation results, we make the above nonlinear partial differential equations (3) in dimensionless form by defining the following parameters: Here, we assume that = = . Equations (3) reduce to the following dimensionless form: where and V represents the dimensionless concentration of substrate and product, and are saturation parameters, and is the reaction diffusion parameter (Thiele modulus). Now, the boundary conditions may be presented as follows [17,18]: = 0, V = 0 when = 0.
The normalized flux becomes

Analytical Expressions for the Concentrations and Current for All Values of Parameters
By using Laplace transform technique and new Homotopy perturbation method (Appendix A), we can obtain the concentrations of substrate and product as follows: where The analytical expression for the dimensionless current is given by where and are defined as in (14).

Analytical Expressions for the Concentrations and Current for Unsaturated (First Order) Kinetics
Now, we consider the limiting case where the substrate concentration is relatively low. In this case, ≤ 1 (i.e., [ ] em ≤ ). Then, (7) will be reduced to the following dimensionless form: The dimensionless form of concentrations obtained by Morf et al. [17] using method of expressions in partial fractions is as follows: where The analytical expression for the dimensionless current is given by (20)

Analytical Solutions for the Concentrations and the Current for Saturated (Zero Order) Kinetics
Next, we consider the limiting case where the substrate concentration is relatively high. In this case, ≥ 1 ([ ] em ≥ ). Equations (7) will be reduced to the following form: The analytical expressions for the concentrations of substrate and product are as follows [22]: The analytical expression for the dimensionless current is given by

Numerical Simulation
The diffusion equations (7) for the corresponding boundary conditions (8), (9), and (10) are solved by numerical methods. The function pdex 4 in Matlab software, which is a function of solving the initial boundary value problems for partial differential equations, was used to solve these equations numerically (Appendix A). The numerical solutions are compared with our analytical results as shown in Figures 1, 2, 3, and 5 and this comparison gives a satisfactory agreement for some possible values of the reaction diffusion parameters.

Results and Discussions
Equations (12) and (13)  of parameters and . The previously reported analytical results ( (17) and (19)) are in terms of the parameter only. Figure 1 shows the time-dependent evolution of normalized concentration profiles for the substrate in the enzyme membrane of a potentiometric sensor. Figures 1(a)-1(c) show dimensionless concentration versus the dimensionless distance . The reaction diffusion parameter is an indicator of the competition between the reaction and diffusion. When is small, the kinetics dominate and the uptake of the substrate are kinetically controlled. From Figure 1(a), it is evident that the value of the substrate concentration decreases when the reaction diffusion parameter increases for different values of . Figure 1(b) illustrates that, when increases, the concentration of the substrate decreases even though the value of is increased. It is obvious from Figure 1(c) that when increases the concentration decreases and if is very small, the concentration of the substrate is uniform and the curve becomes a straight line. Recently, Sivasankari and Rajendran [23] discussed the same mathematical model of potentiometric biosensors for the steady state and according to them when the diffusion parameter is very small, the diffusion of the substrate concentration will be uniform and the curve becomes straight line. This reveals that the parameter time has greater impact on diffusion.
The normalized concentration of the product V for various values of and is plotted in Figures 2(a)-2(c). From the figures we can conclude that the normalized product V increases with the decrease in the value of and increases with the increase in the value of . Figures 3(a) Figure 5(b) illustrates that the dimensionless concentration V also increases with the increase in the dimensionless time . Figure 6(a) characterizes the dimensionless flux versus the dimensionless time for various values of the reaction diffusion parameter and for some and the figure reveals that the value of the flux increases as increases, whereas Figure 6(b) exhibits that the value of flux increases as the value of the parameter decreases.

Figure 5(a) exhibits that the dimensionless concentration increases as the dimensionless time increases and
Our analytical results (12) and numerical results for substrate concentration are compared with previous analytical results (17) of Morf et al. [17] in Figure 7. It exhibits that when is small there is a coincidence between both the results, whereas, when is 1, there is a significant difference between both the results. The same observation for the product concentration is exhibited in Figure 8.

Conclusion
The theoretical analysis of behaviour of potentiometric biosensor was done. The coupled time-dependent nonsteady state nonlinear diffusion equations have been solved analytically and numerically. Moreover we have obtained analytical expressions for the substrate, product concentrations and steady state flux. A good agreement with numerical simulation data is noticed. These analytical results will be used in determining the kinetic characteristics of the biosensor. The theoretical model presented here can be used for the optimization of the design of the biosensor. Furthermore, and where "orange straight line" represents numerical solution, "black dotted line" represents analytical solution of our work (13), and "purple dashed line" represents the analytical solution of Worf 's work (19).
based on the outcome of this work, there is a possibility of extending the procedure to find the approximate amounts of substrate and product concentrations and current for the reciprocal competitive inhibition process.

A. Solution of (7) Using Complex Inversion Formula
In this appendix, we indicate how (12) and (13) are derived.
Using new homotopy perturbation approach [24,25], (7) can be written as (A.1) By using (9) in (A.1), we get where is defined as in (14). Now, by applying Laplace transform to (A.2) and to the conditions in (8), (9), and (10), we obtained the solution of (A.2) as In this appendix, we indicate how (A.4) may be inverted using the complex inversion formula. If ( ) represents the Laplace transform of a function ( ), then, according to the complex inversion formula, we can state that where the integration in (A.5) is to be performed along a line = in the complex plane where = + . The real number is chosen such that = lies to the right of all the singularities but is otherwise assumed to be arbitrary. In practice, the integral is evaluated by considering the contour integral presented on the right-hand side of (A.5), which is then evaluated using the so-called Bromwich contour. The contour integral is then evaluated using the residue theorem which states, for any analytic function ( ), that Hence, we note that (A.11) The first residue in (A.11) is given by The second residue in (A.11) is given by where is defined as in (14). Here, we used cosh( ) = cos( ) and sinh( ) = sin( ). From (A.11), (A.12), and (A.13), we conclude that where is defined as in (14). Similarly, we can solve (A.3) by using complex inversion formula.