#### Abstract

Mathematical modeling of amperometric biosensor with cyclic reaction is discussed. Analytical expressions pertaining to the concentration of substrate, cosubstrate, reducing agent and medial product and current for hybrid enzyme biosensor are obtained in terms of Thiele module and saturation parameters. In this paper, a powerful analytical method, called homotopy analysis method (HAM) is used to solve the system of nonlinear differential equations. Furthermore, in this work the numerical simulation of the problem is also reported using Scilab/Matlab program. Our analytical results are compared with simulation results. A good agreement between analytical and numerical results is noted.

#### 1. Introduction

Biosensor (Figure 1) is a device that uses specific biochemical reactions mediated by isolated enzymes, immunosystems, tissues, organelles, or whole cells to detect chemical compounds usually by electrical, thermal, or optical signals [1]. They involve a biological (recognition) element and a transduction element. The biological or recognition element may be an antibody, an enzyme, DNA, RNA, a whole cell, or a whole organ or system. The transduction element, wherein the biological event or signal is converted to a measurable signal, may include anyone of the following forms: chemical, electrical, magnetic, mechanical, optical, or thermal.

The biosensor was first described by Clark and Lyons in 1962, when the term enzyme electrode was adopted [2]. The term “biosensor” was introduced by Cammann in 1977 [3]. Since then, research communities from various fields such as physics, chemistry, and material science have come together to develop more sophisticated, reliable, and mature biosensing devices for applications in the fields of medicine, agriculture, biotechnology, as well as in the military for bioterrorism detection and prevention [4]. Biosensors offer the prospects of simplified, virtually nondestructive analysis of turbid biological fluids. Also, biosensors for medical care have demanded the greatest attention for technical development [5].

Amperometric electrodes have been used in the design of biosensors for glucose, aminoacids, and other molecules [6–9]. In cases of amperometric enzyme biosensors the potential at the electrode is held constant while the current flow is measured. Amperometric biosensors are quite sensitive and more suited for mass production than the potentiometric ones [10, 11]. Electropolymerized films offer wide immobilization capabilities and extremely large diversity in the development of biosensors [12]. The development of amperometric biosensors is an area of growing interest in many branches of science [13–17]. Hybrid biosensors are biosensors with more than one biosensitive material-enzyme, tissue microorganism [18] or other agents [19]. Hybrid biosensors are less selective than enzymatic ones, but the choice of the type of biosensor should be made according to several additional factors like availability of biocatalysts, their stability, and sensitivity required [20].

Recently, Rajendran and coworkers [21–23] obtained the analytical expressions of substrate and product for amperometric biosensors. Eswari and Rajendran [24, 25] obtained the steady-state current at microdisk biosensor and microcylinder biosensor. Recently, Rangelova [26] found the concentration profiles of substrate, cosubstrate, reducing agent, medial product using finite difference technique and Matlab program. To our knowledge, no general simple analytical expressions for the concentrations and current for hybrid biosensors have been reported for all values of the Thiele module and saturation parameters. However, in general, analytical solution of nonlinear differential equations is more interesting and useful than purely numerical solutions as they are amenable to various kinds of manipulation and data analysis. Analytical expressions are usually derived from the basic physical principles and free from numerical dispersions and other truncation errors that often occurred in numerical simulations. For this reason, in this paper we have derived an analytical expression for the concentrations and current for hybrid enzyme biosensor for all values of the normalized parameters using homotopy analysis method (HAM).

#### 2. Mathematical Formulation and Analysis of the Problem

##### 2.1. Mathematical Formulation

In enzyme-based catechol biosensor, the cyclic reaction scheme for the substrate, co-substrate, reducing agent and medial product can be represented as follows [26]: where and are reaction rate constants, is the first product and is the second product (dehydroascorbic acid), is the measured substrate (catechol), is the cosubstrate (oxygen), is the reducing agent (-ascorbic acid), and is the medial product (1, 2 benzoquinone).

A cyclic reaction between catechol and 1,2-benzoquinone takes place by combining the tyrosinase reaction and the chemical reduction of 1,2 benzoquinone to catechol by -ascorbic acid. This acid is well known as an effective reducing agent, and 1,2 benzoquinone would be reduced to catechol and may drive the reaction in the opposite direction to that of enzymatic oxidation of catechol. The oxygen consumed in the enzymatic reaction is not compensated for by the chemical reduction. Therefore, if -ascorbic acid does not affect the enzyme activity of tyrosinase, a cyclic reaction should take place and the consumption of the dissolved oxygen will continue until its concentration becomes zero [27]. The differential equations for this ping-pong mechanism at the steady-state condition are as follows [26]:

Here (mM) denote the concentration of substrate, cosubstrate, reducing agent and medial product, respectively. (m^{2}/s) are diffusion coefficients, (mmol/(l.s)) are the reaction rate constants, is the distance coordinate, (mM) are the reaction rate constants, and (mmol/(l.s)) is the maximal rate. The boundary conditions are given by [26]
where is the active membrane thickness. The diffusion limiting current can be expressed as follows:
where is the number of electrons taking part in the electrochemical reaction on the cathode, is the Faraday’s constant, and is the surface of cathode.

##### 2.2. Normalized Form

By introducing the following set of nondimensional variables, the nonlinear reaction/diffusion equation (2) take the following normalized form: The transformed boundary conditions are The dimensionless form of the current is given by

#### 3. Approximate Analytical Expression of Concentrations of Four Reactants Using Homotopy Analysis Method

##### 3.1. Homotopy Analysis Method

The homotopy analysis method (HAM) [28–36] is a general analytic method to get series solutions of various types of nonlinear equations, including ordinary differential equations, partial differential equations, and coupled nonlinear equations. Unlike perturbation methods, the HAM is independent of small/large physical parameters. More importantly, different from all perturbation and traditional nonperturbation methods, the HAM provides us with a simple way to ensure the convergence of solution series. Besides, different from all perturbation and previous nonperturbation methods, the HAM provides us with great freedom to choose proper base functions to approximate a nonlinear problem [29, 34]. Liao’s book [29] for the homotopy analysis method was first published in 2003. Now, more and more researchers have been successfully applying this method to various nonlinear problems in science and engineering. In this paper we employ HAM to solve the nonlinear differential equations (6), (7), and (9). The basic concept of Homotopy analysis method is given in Appendix A.

##### 3.2. Solution of Boundary Value Problem Using the Homotopy Analysis Method

Using HAM method (Appendix B), we obtain the analytical expression corresponding to the concentrations of substrate, co-substrate, reducing agent and medial product as follows: where , , and is the convergence control parameter. Equation (12)–(15) represent the new analytical expression of the dimensionless reactant concentrations. Using (11) and (15), we can obtain the current as follows:

#### 4. Result and Discussion

##### 4.1. Numerical Simulation

In order to investigate the accuracy of this analytical method with a finite number of terms, the system of differential equations ((6)–(9)) also solved by numerical methods. The function pdepe (finite element method) in Scilab/Matlab software which is a function of solving PDE is used to solve these nonlinear equations [37]. The Scilab/Matlab program is also given in Appendix C. To validate the results, the convergence studies are carried out. The convergence region of auxiliary parameter is given in the Appendix D. To show the efficiency of the present analytical method, our results are compared with the numerical solution (Scilab/Matlab program) in Tables 1, 2, 3, and 4 and Figures 2–6. The average relative errors between our analytical results and numerical results are 0.23%, 0.11%, 0.03%, and 0.01% for the concentrations of substrate, co-substrate, reducing agent and medial product, respectively.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

##### 4.2. Effect of the Thiele Module

The concentration of substrate, co-substrate, reducing agent and medial product depends upon Thiele module and saturation parameters. The Thiele module essentially compares the rate of enzyme reaction (/) and diffusion in the enzyme layer (). We observe the rise and downfall of concentration profiles in two cases. (i) If Thiele module is small (), then enzyme kinetics predominates in the biosensor response. The overall kinetics is governed by the total amount of active enzyme. (ii) The response is under diffusion control, if the Thiele module is large , which is observed at high catalytic activity and active membrane thickness or at low reaction rate constant or diffusion coefficient values .

The concentration profiles for the four reactants for some fixed values of parameters are shown in Figure 2. The analytical results are compared with the numerical result, graphically. Upon comparison, it is evident that both results give satisfactory agreement. The concentration of substrate increases and attains its maximum, when Thiele module and saturation parameter increase (refer to Figure 3). Therefore the profile deviates more from the linearity. The concentration profile representing medial product increases as Thiele module and saturation parameter increase (refer to Figure 4). From Figure 5, it is inferred that the increasing value of coordinate decreases the concentration of reducing agent . Figure 6 represents the concentration of co-substrate versus the normalized distance coordinate. It reaches the minimum value zero at electrode interface.

##### 4.3. Influence of Active Membrane Thickness

The active membrane thickness is one of the important technical parameters and it has a considerable effect on the concentration profiles for the four reactants. Also, the Thiele module is directly proportional to the thickness of the membrane. The concentration profiles depend significantly on the membrane thickness . If the active membrane thickness is large ( or thick membrane), the concentration profiles of substrate and medial product increase whereas the concentration of reducing agent and co-substrate decreases. If the membrane thickness is small ( or thinner membrane), the concentration profiles of substrate, medial product and reducing agent have uniform values. These uniform values are equal to the concentration of the above reactants at . The concentration profile of co-substrate increases when thickness of the membrane decreases. The increasing value is not significant. Also, the concentration of co-substrate increases when the saturation parameter increases (refer to Figures 3–6).

##### 4.4. Current Response

The normalized current versus the diffusion coefficient ratio is calculated at different values of the active membrane thickness . The results obtained for various values of the normalized parameters are depicted in Figure 7(a). The current response increases when the active membrane thickness () increases. Also, for thinner membrane () the value of the current is zero. In Figure 7(b), the current response increases, as the saturation parameter increases.

**(a)**

**(b)**

#### 5. Conclusions

The theoretical model of hybrid amperometric enzyme biosensor with cyclic reaction and biochemical amplification for steady-state condition is discussed. The system of three nonlinear differential equations for ping-pong enzyme kinetics has been solved analytically. Influence of Thiele module and active membrane thickness is investigated. The obtained results have a good agreement with those obtained using numerical method. This analytical result will be useful in sensor design, optimization, and prediction of the electrode response. Using this result, the action of biosensor is analyzed at critical concentration of substrate and enzyme activities. Theoretical results obtained in this paper can also be used to analyze the effect of different parameters such as active membrane thickness and saturation parameters.

#### Appendices

#### A. Basic Idea of Liao’s Homotopy Analysis Method

Consider the following differential equation [38]:
where is a nonlinear operator, denotes an independent variable, and is an unknown function. For simplicity, we ignore all boundary or initial conditions, which can be treated in a similar way. By means of generalizing the conventional Homotopy method, Liao constructed the so-called zero-order deformation equation as
where is the embedding parameter, is a nonzero auxiliary parameter, is an auxiliary function, is an auxiliary linear operator, is an initial guess of , and is an unknown function. It is important that one has great freedom to choose auxiliary unknowns in HAM. Obviously, when and, it holds
respectively. Thus, as increases from 0 to 1, the solution varies from the initial guess to the solution . Expanding in Taylor series with respect to *,* we have
where
If the auxiliary linear operator, the initial guess, the auxiliary parameter *,* and the auxiliary function are so properly chosen, the series (A.4) converges at then we have
Define the vector
Differentiating (A.2) for times with respect to the embedding parameter , then setting , and finally dividing them by !, we will have the so-called th-order deformation equation as
where
Applying on both sides of (A.8), we get
In this way, it is easy to obtain for , at th order; we have
When , we get an accurate approximation of the original equation (A.1). For the convergence of the previous method we refer the reader to Liao [29]. If (A.1) admits unique solution, then this method will produce the unique solution. If (A.1) does not possess unique solution, the HAM will give a solution among many other (possible) solutions.

#### B. Approximate Analytical Expression of Concentrations of Substrate, Co-Substrate, Reducing Agent and Medial Product

From (8) it is clear that the concentration of normalized reducing agent is In order to solve (6), (7), and (9) by means of the HAM, we first construct the zeroth-order deformation equation by taking : The approximate solutions of (B.2)–(B.4) are as follows: Substituting (B.5) in (B.2), (B.6) in (B.3), and (B.7) in (B.4) and equating the like powers of we get The boundary conditions equation (10) become From (B.8), (B.10), and (B.12) and from the boundary conditions (B.14) and (B.15), we get Substituting the values of in (B.9), (B.11), and (B.13) and solving the equations using the boundary conditions (B.16) and (B.17) we obtain the following results: where

Adding (B.18) and (B.21), we get (12) in the text. Similarly we get (13) and (15) in the text.

#### C. Scilab/Matlab Program to Find the Numerical Solution of Nonlinear Equations (6)–(9)

function pdex4 m = 0; x = linspace(0,1); t=linspace(0,100000); sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t); u1 = sol(:,:,1); u2 = sol(:,:,2); u3 = sol(:,:,3); u4 = sol(:,:,4); figure plot(x,u1(end,:)) title(‘u1(x,t)’) xlabel(‘Distance x’) ylabel(‘u1(x,2)’) figure plot(x,u2(end,:)) title(‘u2(x,t)’) xlabel(‘Distance x’) ylabel(‘u2(x,2)’) figure plot(x,u3(end,:)) title(‘u3(x,t)’) xlabel(‘Distance x’) ylabel(‘u3(x,2)’) figure plot(x,u4(end,:)) title(‘u4(x,t)’) xlabel(‘Distance x’) ylabel(‘u4(x,2)’) function [c,f,s] = pdex4pde(x,t,u,DuDx) c = [1; 1; 1; 1]; f = [1; 1; 1; 1].* DuDx; Q=9; p1=0.92; p=1.2; m1=0.2; m2=0.05; m3=0.001; n1=0.8; n2=0.1; n3=0.6; F=-Q*p1*((1/(1+1/(u(1))+1/(u(4)))-m1*(u(2))* (u(3)))); F1=-Q*n2*(m2*(u(2))-1/(1+1/(u(1))+1/(u(4)))); F2=-m3*n3*Q*(u(3)); F3=-Q*n1*p*(1/(1+1/(u(1))+1/(u(4)))); s=[F; F1; F2; F3]; function u0 = pdex4ic(x); u0 = [1; 1; 1; 1]; function [pl,ql,pr,qr]=pdex4bc(xl,ul,xr,ur,t) pl = [ul(1)-0.09;ul(2)−0;ul(3)-1;ul(4)-0.54]; ql = [0; 0; 0; 0]; pr = [0; 0; 0; ur(4)]; qr = [1; 1; 1; 0];

#### D. Determining the Validity Region of

The analytical solution represented by (12), (13), and (15) contains the auxiliary parameter , which gives the convergence region and rate of approximation for homotopy analysis method. The analytical solution should converge. It should be noted that the auxiliary parameter controls the convergence and accuracy of the solution series. In order to define region such that the solution series is independent of , a multiple of curves are plotted. The region where the distribution of versus is a horizontal line is known as the convergence region for the corresponding function. The common region among is known as the overall convergence region. To study the influence of on the convergence of solution, the curves of are plotted in Figure 8. This figure clearly indicates that the valid region of is about . Similarly we can find the value of the convergence-control parameter for different values of constant parameters.

#### Nomenclature

*Symbols*

: | Measured substrate concentration of catechol (mM) |

: | Medial product concentration of 1,2 benzoquinone (mM) |

: | Reducing agent concentration of L-ascorbic acid (mM) |

: | Co-substrate concentration of oxygen (mM) |

: | Maximal rate (mmol/(l.s)) |

: | Diffusion coefficient for substrate (/s) |

: | Diffusion coefficient for medial product (/s) |

: | Diffusion coefficient for reducing agent (/s) |

: | Diffusion coefficient for co-substrate (/s) |

: | Reaction rate constants (mmol/(l.s)) |

: | Reaction rate constants (mM) |

: | Distance coordinate () |

: | Active membrane thickness () |

: | Convergence control parameter |

: | Normalized measured substrate concentration (dimensionless) |

: | Normalized medial product concentration (dimensionless) |

: | Normalized reducing agent concentration (dimensionless) |

: | Normalized co-substrate concentration (dimensionless) |

: | Normalized distance coordinate (dimensionless) |

: | Saturation parameters (dimensionless) |

: | Linear enzyme kinetic coefficient (dimensionless) |

: | Ratio of diffusion coefficients (dimensionless) |

: | Ratio of reaction rate constants (dimensionless) |

: | Thiele module (dimensionless) |

: | Normalized current (dimensionless). |

#### Acknowledgments

This work was supported by the University Grants Commission (F. no. 39–58/2010(SR)), New Delhi, India. The authors are thankful to Dr. R. Murali, The Principal, The Madura College, Madurai, and Mr. M. S. Meenakshisundaram, The Secretary, Madura College Board, Madurai, for their encouragement. The author K. Indira is very thankful to the Manonmaniam Sundaranar University, Tirunelveli, for allowing to do the research work.