Theoretical Analysis of an Amperometric Biosensor Based on Parallel Substrates Conversion

Abstract

The behaviour of an amperometric biosensor based on parallel substrates conversion for steady-state condition has been discussed. This analysis contains a nonlinear term related to enzyme kinetics. Simple and closed form of analytical expressions of concentrations and of biosensor current is derived. This model was originally reported by Vytautas Aseris and his team (2012). Concentrations of substrate and product are expressed in terms of single dimensionless parameter. A new approach to Homotopy perturbation method (HPM) is employed to solve the system of nonlinear reaction diffusion equations. Furthermore, in this work, the numerical solution of the problem is also reported using Matlab program. The analytical results are compared with the numerical results. The analytical result provided is reliable and efficient to understand the behavior of the system.

1. Introduction

A biosensor is an analytical device, used for the detection of an analyte, which combines a biological component with a physiochemical detector [1]. Analyzed substrate is biochemically converted to a product when the biosensor operation is being undertaken. In most of the cases, biosensor response is directly proportional to the concentration of the reaction product [2]. The amperometric biosensors have proved to be reliable in various systems with applications in the field of medicine, food technology, and the environmental industry [3, 4]. The understanding of the kinetic peculiarities of biosensors is of crucial importance for their design. The mathematical modeling is rather widely used to improve the efficiency of the biosensors design and to optimize their configuration [5–8]. Since 1970s, various mathematical models of biosensors have been developed and used successfully to study and optimize analytical characteristics of biosensors [9–14]. Mathematical modeling of two-enzyme biosensors has been started in 1980s with the modeling of an amperometric monolayer enzyme electrode with two coimmobilized enzymes [15, 16]. Later, nonlinear mathematical models have been developed for amperometric two-enzyme biosensors with different enzymes [17, 18].

The numerical method of solving partial differential equations is to make calculation at all intervals of substrates concentration and at different diffusion and enzymatic reaction rates. The nonstationary diffusion equations [19], containing a nonlinear term related to the enzymatic reaction, are carried out using the implicit difference scheme [20]. In recent years, analytical solutions are reported for various types of biosensors [21–27].

To the best of our knowledge, there is no simple analytical expression for steady-state concentrations and current of amperometric-based biosensor utilizing parallel substrates conversion for all values of parameters. In this paper, we present the simple and closed form of the analytical expressions for the concentrations of substrate and of product using new Homotopy perturbation method (HPM) [28]. The current density was determined corresponding to all possible values of the parameters.

2. Mathematical Formulation of the Problem and Analysis

The reaction scheme employing two parallel substrates (𝑆₁ and 𝑆₂), two enzymes (𝐸₁ and 𝐸₂), and two reaction products (𝑃₁ and 𝑃₂) in biosensor may be as follows: 𝑆→₁^𝐸112𝑃₁,𝑆₁+𝑆→₂^𝐸2𝑃₂. (1)

First phase of the biosensor action contains hydrogen peroxide (𝑆₁) only. The biosensor response reaches the steady-state at the end of the first phase. The second phase of the biosensor action begins at the second substrate (𝑆₂) poured into the solution. The relative changes between the responses in both phases are measured as the final response of the biosensor. The schematic diagram of the system is provided in Figure 1. The above reaction in the enzymatic layer with the one-dimensional-in-space diffusion, described by Fick’s law, leads to the following nonlinear transient reaction-diffusion equations [29]: 𝐷_𝑆1,𝑒𝜕²𝑆_1,𝑒𝜕𝑥²−𝑘₁𝑒₁𝑆_1,𝑒−𝑘₂₁𝑘₂₂𝑒₂𝑆_1,𝑒𝑆_2,𝑒𝑘₂₁𝑆_1,𝑒+𝑘₂₂𝑆_2,𝑒𝐷=0,_𝑆2,𝑒𝜕²𝑆_2,𝑒𝜕𝑥²−𝑘₂₁𝑘₂₂𝑒₂𝑆_1,𝑒𝑆_2,𝑒𝑘₂₁𝑆_1,𝑒+𝑘₂₂𝑆_2,𝑒𝐷=0,_𝑃1,𝑒𝜕²𝑃_1,𝑒𝜕𝑥²+𝑘₁𝑒₁𝑆_1,𝑒2𝐷=0,_𝑃2,𝑒𝜕²𝑃_2,𝑒𝜕𝑥²+𝑘₂₁𝑘₂₂𝑒₂𝑆_1,𝑒𝑆_2,𝑒𝑘₂₁𝑆_1,𝑒+𝑘₂₂𝑆_2,𝑒=0, (2) where 𝑥 denotes space, 𝑆_1,𝑒, 𝑆_2,𝑒 are the molar concentration of the substrates, and 𝑃_1,𝑒, 𝑃_2,𝑒 are the molar concentration of products. 𝐷_𝑆1,𝑒, 𝐷_𝑆2,𝑒, 𝐷_𝑃1,𝑒, and 𝐷_𝑃2,𝑒 are the corresponding diffusion coefficients and 𝑘₁, 𝑘₂₁, and 𝑘₂₂ are the reaction rate constants. Graphical representation of the products 𝑃_1,𝑒, 𝑃_2,𝑒 and substrates 𝑆_1,𝑒, 𝑆_2,𝑒, with steady-state boundary conditions employed in the simulated multilayered biosensor is provided in Figure 2. The mass transport by diffusion of the substrates and products takes place outside the enzyme layer. The diffusion layer is treated as the Nernst diffusion layer. The boundary conditions for the above system of nonlinear equation are as follows: 𝐷_𝐶𝑒𝜕𝑐_𝑒𝜕𝑥_𝑥=0=0,𝑃_1,𝑒(𝑥=0)=0,𝑐=𝑆₁,𝑆₂,𝑃₂𝑆_𝑖,𝑒𝑥=𝑑_𝑒+𝑑_𝑑=𝑆_𝑖0,𝑃_𝑖,𝑒𝑥=𝑑_𝑒+𝑑_𝑑=0,𝑖=1,2, (3) where 𝑑_𝑒 and 𝑑_𝑑 denote the thickness of the enzyme and diffusion layer. For the steady-state condition, the following sets of dimensionless variables are introduced: 𝑆_1𝑁=𝑆_1,𝑒𝑆₁₀,𝑆_2𝑁=𝑆_2,𝑒𝑆₂₀,𝑃_1𝑁=𝑃_1,𝑒𝑆₁₀,𝑃_2𝑁=𝑃_2,𝑒𝑆₁₀,𝑥𝑦=𝑑,𝛾_𝑆1=𝑑²𝑘₁𝑒₁𝐷_𝑆1,𝑒,𝛾_𝑆2=𝑑²𝑘₂₂𝑒₂𝐷_𝑆1,𝑒,𝛾_𝑆3=𝑑²𝑘₂₂𝑒₂𝐷_𝑆2,𝑒,𝛾_𝑃1=𝑑²𝑘₁𝑒₁𝐷_𝑃1,𝑒,𝛾_𝑃2=𝑑²𝑘₁𝑒₁𝐷_𝑃2,𝑒,𝑆𝛼=¹⁰𝑆₂₀𝑘,𝛽=²²𝑘₂₁, (4) where 𝑑=𝑑_𝑒+𝑑_𝑑. The governing nonlinear reaction diffusion equations (2) are expressed in the following steady-state nondimensional form: 𝑑²𝑆_1𝑁𝑑𝑦²=𝛾_𝑆1𝑆_1𝑁+𝛾_𝑆2𝑆_1𝑁𝑆_2𝑁𝛼𝑆_1𝑁+𝛽𝑆_2𝑁,𝑑²𝑆_2𝑁𝑑𝑦²=𝛾_𝑆3𝑆_1𝑁𝑆_2𝑁𝛼𝑆_1𝑁+𝛽𝑆_2𝑁,𝑑²𝑃_1𝑁𝑑𝑦²=−𝛾_𝑃1𝑆_1𝑁2,𝑑²𝑃_2𝑁𝑑𝑦²=−𝛾_𝑃2𝑆_1𝑁𝑆_2𝑁𝛼𝑆_1𝑁+𝛽𝑆_2𝑁, (5)(6)(7)(8) where 𝑆_1𝑁 and 𝑆_2𝑁 represent the dimensionless substrate concentrations, 𝑃_1𝑁 and 𝑃_2𝑁 represent the dimensionless product concentrations, and 𝑦 represents the dimensionless distance. The reaction/diffusion parameters of substrates and products are considered to be 𝛾_𝑆1, 𝛾_𝑆2, 𝛾_𝑆3, 𝛾_𝑃1, and 𝛾_𝑃2. The boundary conditions are represented as follows: 𝑑𝑆_1𝑁=𝑑𝑦𝑑𝑆_2𝑁=𝑑𝑦𝑑𝑃_2𝑁𝑑𝑦=0,𝑃_1𝑁𝑆=0when𝑦=0,_1𝑁=𝑆_2𝑁=1,𝑃_1𝑁=𝑃_2𝑁=0when𝑦=1. (9) The dimensionless current density is given by 𝑖𝜓=^𝑠𝑑𝑛_𝑒𝐹𝐷_𝑃1,𝑒𝐴𝑆₁₀=𝑑𝑃_1𝑁𝑑𝑦_𝑦=0. (10)

Figure 1 

Schematic representation of biosensor utilising parallel substrate conversion with enzymatic layer (of thickness 𝑑𝑒) and the diffusional layer (of thickness 𝑑𝑑) the bioelectrode functioning in the presence of hydrogen peroxide (𝑆1) or catalyses (𝑆2). 𝐸1 and 𝐸2 represent the catalase and peroxidase activities. 𝑃1 and 𝑃2 are products.

Figure 2 

Graphical representation of steady-state boundary conditions of the products 𝑃1,𝑒, 𝑃2,𝑒 and substrates 𝑆1,𝑒, 𝑆2,𝑒, for multilayered biosensor. Here, 𝑥=0 is taken at the surface of the electrode and 𝑥=𝑑𝑒+𝑑𝑑 is taken as the surface of the diffusional layer.

3. A New Approach to Homotopy Perturbation Method and Its Advantages

The Homotopy perturbation method is first introduced to solve nonlinear wave equations. The advantage of the method is that it does not need a small parameter in the system, leading to wide application in nonlinear wave equations [30]. Recently, many authors have used HPM to various problems and reported the efficiency of the HPM for handling nonlinear engineering problems [31–34]. This method is a combination of Homotopy in topology and classic perturbation techniques. He used the HPM to solve the Light hill equation [35]; this wide variety of applications shows the power of the HPM to solve functional equations. The HPM has uniqueness in its applicability, accuracy, and efficiency. Recently, a new approach to HPM is introduced to solve the nonlinear problem, in which one will get better simple approximate solution in the zeroth iteration [28].

4. Analytical Expression for the Concentration of Substrate and Product Using New Homotopy Perturbation Method

In this paper, a new approach to Homotopy perturbation method is used (Appendix A) to solve the nonlinear differential equations (5)–(8). Using this method, the analytical expression of the two substrate concentrations (𝑆_1𝑁 and 𝑆_2𝑁) and two product concentrations (𝑃_1𝑁 and 𝑃_2𝑁) can be obtained as follows:𝑆_1𝑁𝑦=cosh𝐴𝑦(𝐴),cosh𝑆_2𝑁𝑦=cosh𝐵𝑦(𝐵),cosh𝑃_1𝑁𝑦=𝛾_𝑃11−cosh𝐴𝑦+𝑦cosh(𝐴)−𝑦2𝐴²cosh(𝐴)𝑃_2𝑁𝑦=𝛾_𝑃2𝛾_𝑆31−cosh𝐵𝑦(𝐵),cosh (11)(12)(13)(14) where 𝐴=𝛾_𝑆1𝛼+𝛽+𝛾_𝑆2𝛼+𝛽,𝐵=𝛾_𝑆3.𝛼+𝛽 (15) Equations (11)–(14) satisfy the boundary conditions (9). These equations represent the simple and closed form of analytical expressions of concentrations of substrate and product for all possible values of the parameters. Using (6) and (8) with boundary conditions (9), the following relation is also obtained: 𝑃_2𝑁𝑦=𝛾_𝑃2𝛾_𝑆31−𝑆_2𝑁𝑦. (16) The analytical expression of the current using (10) is as follows: 𝜓=𝑖𝑑𝑛_𝑒𝐹𝐴𝐷_𝑃1,𝑒𝑆₁₀=𝛾_𝑃1[cosh(𝐴)−1]2𝐴²(𝐴).cosh (17) When 𝐴→0, the above expression can be simplified as 𝜓=𝛾_𝑃1/4.

4.1. Limiting Case

The mathematical as well as the corresponding computational models of the biosensor were evaluated. Assume that zero concentration of the second substrate (𝑆₂₀=0) leads to reducing the governing equations (5)–(8) to the following equations: 𝑑²𝑆_1𝑁𝑑𝑦²=𝛾_𝑆1𝑆_1𝑁,𝑑²𝑃_1𝑁𝑑𝑦²=−𝛾_𝑃1𝑆_1𝑁2. (18) The steady-state solution of the above equations for the boundary condition (9) is given by 𝑆_1,𝑁=cosh√𝛾_𝑆1𝑦cosh√𝛾_𝑆1,𝑃_1𝑁𝑦=𝛾_𝑃11−cosh√𝛾_𝑆1𝑦+𝑦cosh√𝛾_𝑆1−𝑦2𝛾_𝑆1cosh√𝛾_𝑆1. (19) The dimensionless steady-state current using (10) is given by 𝜓=𝑖𝑑𝑛_𝑒𝐹𝐴𝐷_𝑃1,𝑒𝑆₁₀=𝛾_𝑃1cosh√𝛾_𝑆1−12𝐴²cosh√𝛾_𝑆1. (20) Equation (20) is identical to (17), when 𝑆₂₀=0 (i.e., 𝛼=∞) in (17).

5. Comparison of Analytical Expressions with Experimental Work and Numerical Simulation

To show the efficiency, the steady-state result (11)–(14) of the reaction diffusion equations is compared with the numerical solution of (5)–(8). The numerical solutions are found using pdepe (finite element method) in Matlab [36]. This Matlab program [37] is also given in Appendix B. The comparison is shown in Figures 3–6 and it gives satisfactory agreement.

Figure 3 

Normalized substrates and products concentration versus dimensionless space 𝑦 were computed using (11)–(14) for various values of reaction/diffusion parameters and for some fixed values of the other existing parameters. The key to the plot: (—) (11)–(14) and (∙ ∙∙) and (+++) numerical simulation of (5)–(8).

(a)

(b)

(c)

(d)

Figure 4 

Normalized substrate concentration 𝑆1𝑁 versus dimensionless space 𝑦 was computed using (11) for various values of reaction/diffusion parameter 𝛾𝑆1, dimensionless parameter 𝛼 and for some fixed values of the parameters 𝛾𝑆1, 𝛾𝑆2, 𝛼, 𝛽. The key to the plot: (—) (11) and (∙ ∙∙) numerical simulation of (5).

(a)

(b)

(c)

(d)

Figure 5 

Normalized substrate concentration 𝑆2𝑁 versus dimensionless space 𝑦 was computed using (12) for various values of reaction/diffusion parameter 𝛾𝑆3, dimensionless parameter 𝛼 and for some fixed values of the parameters 𝛼, 𝛽, 𝛾𝑆3. The key to the plot: (—) (12) and (∙ ∙∙) numerical simulation of (6).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6 

Normalized product concentration 𝑃1𝑁 versus dimensionless space 𝑦 was computed using (13) for various values of reaction/diffusion parameter 𝛾𝑃1, 𝛾𝑆1 and for some fixed values of the parameters 𝛾𝑆1, 𝛾𝑆2, 𝛾𝑃1, 𝛼, 𝛽. The key to the plot: (—) (13) and (∙ ∙∙) numerical simulation of (7).

6. Discussion

Equations (11)–(14) are the new simple and closed approximate analytical expressions for the dimensionless substrate concentrations (𝑆_1𝑁 and 𝑆_2𝑁) and product concentrations (𝑃_1𝑁 and 𝑃_2𝑁) in terms of 𝐴 or 𝐵. The analytical results (11)–(14) are compared with numerical solution, which is provided in Figures 3–7. Figure 3 shows the plots of all the concentration versus dimensionless distance 𝑦.

(a)

(b)

(c)

(d)

Figure 7 

Normalized product concentration 𝑃2𝑁 versus dimensionless space 𝑦 was computed using (14) for various values of reaction/diffusion parameter 𝛾𝑃2, dimensionless parameter 𝛼 and for some fixed values of the parameters 𝛾𝑃1, 𝛾𝑆3, 𝛼, 𝛽. The key to the plot: (—) (14) and (∙ ∙∙) numerical simulation of (8).

The kinetic response of a biosensor depends on the concentration species. The concentration depends on the following factors: 𝛾_𝑆1, 𝛾_𝑆2, 𝛾_𝑆3, 𝛾_𝑃1, and 𝛾_𝑃2. The diffusion parameters 𝛾_𝑆1 and 𝛾_𝑆3 represent the ratio of the characteristic time of the enzymatic reaction to that of concentration diffusion. These parameters can be varied by changing the thickness of the enzyme or diffusional layer. These parameters describe the relative importance of diffusion and reaction in both layers. When 𝛾_𝑆1 or 𝛾_𝑆3 is small, the kinetics is dominant resistance; the uptake of concentration species is kinetically controlled. Under these conditions, the concentration species across both the layers is essentially uniform. When the diffusion parameter 𝛾_𝑆1 or 𝛾_𝑆3 is large, diffusion limitations are the principal determining factor.

Figures 4 and 5 show the dimensionless concentration of substrates (𝑆_1𝑁 and 𝑆_2𝑁) versus dimensionless distance 𝑦. Figures 3 and 4(a) and 4(b) are plotted for various values of diffusion parameters 𝛾_𝑆1 and 𝛾_𝑆3 for some fixed values of 𝛼,𝛽 and 𝛾_𝑆2. From these figures, it is clear that the concentration of substrates 𝑆_1𝑁 and 𝑆_2𝑁 decreases when the value of 𝛾_𝑆1 and 𝛾_𝑆3 increases, respectively. The concentration is approximately equal to one when 𝛾_𝑆1≤0.01 and 𝛾_𝑆3≤0.01. Figures 4(c) and 5(c) are plotted for various values of dimensionless parameter 𝛼 for some fixed values of other existing parameters. From these figures, it is clear that concentration increases for large values of 𝛼. Figures 4(d) and 5(d) are plotted for various values of 𝐴 and 𝐵, respectively. It is evident that the concentration is uniform, when 𝛾_𝑆1(𝛼+𝛽)+𝛾_𝑆2<(0.01)(𝛼+𝛽) (or 𝐴<0.1) and 𝛾_𝑆3<(0.01)(𝛼+𝛽) (or 𝐵<0.1).

Figures 6(a)–6(f) show the dimensionless product concentration profile 𝑃_1𝑁 versus dimensionless distance 𝑦 using (13). Figures 6(a) and 6(b) are plotted for all various values of 𝛾_𝑃1 for some fixed values of 𝛼, 𝛽, 𝛾_𝑆1 and 𝛾_𝑆2. Thus, it is concluded that there is a simultaneous increase in the values of the concentration 𝑃_1𝑁 as well as in 𝛾_𝑃1. Figure 6(c) is plotted for various values of 𝛼 and the fixed value of 𝛾_𝑆1, 𝛾_𝑆2, 𝛾_𝑃1, and 𝛽. It shows that the concentration is always minimum for all values of 𝛼. Concentration versus distance is plotted in Figures 6(d) and 6(e) for various values of 𝛾_𝑆1 and for some fixed values of other existing parameters. Thus, it is concluded that the concentration is always minimum and also there is a simultaneous decrease in the values of the concentration 𝑃_1𝑁 as the value of 𝛾_𝑆1 increases. Figure 6(f) is plotted for various values of 𝐴, and from these graphs, it is inferred that concentration is minimum and uniform when 𝐴≥0.

Figures 7(a)–7(d) are plotted for the dimensionless product concentration 𝑃_2𝑁 as ordinate versus dimensionless space 𝑦 as abscissa using (14). From Figures 7(a) and 7(b), it is evident that the concentration of 𝑃_2𝑁 increases as the diffusivity parameter 𝛾_𝑃2 increases. The concentration is approximately equal to zero, that is, (𝑃_2𝑁≈0) when 𝛾_𝑃2≤0.1 for all possible values of the parameters 𝛼, 𝛽 and 𝛾_𝑆3. From Figure 7(c), it is inferred that concentration 𝑃_2𝑁 decreases, when the value of 𝛼 increases and 𝑃_2𝑁≈0 when 𝛼≥10. Figure 7(d) is plotted for various values of 𝐵; it is evident that the concentration increases as the value of 𝐵 increases. The concentration increases, when 𝛾_𝑆3>(0.01)(𝛼+𝛽) or 𝐵>0.1.

Figures 8(a)–8(e) are plotted using (17), which shows the dynamics of normalized steady-state current 𝜓 versus dimensionless parameters 𝛾_𝑃1, (𝛼+𝛽), 𝛾_𝑆1, 𝛾_𝑆2, and 𝐴, respectively. From Figure 8(a), it is noted that, for small values of 𝛾_𝑆1(𝛾_𝑆1≤0.1) and for any values of 𝛼 and 𝛽, the value of the current is significantly high. The current grows notably faster at higher values of 𝛾_𝑆1 but in the case of 𝛾_𝑆1≥10 and 𝛼=𝛽=0.1, the current does not increase. Figure 8(b) is plotted for various values of 𝛾_𝑃1 and for some fixed values of 𝛾_𝑆1 and 𝛾_𝑆2. It is observed that the current increases rapidly and reaches the steady-state when (𝛼+𝛽)≥0.5. From Figures 8(c) and 8(d), it is evident that the current is uniformly decreasing for large values of 𝛾_𝑃1, that is, (𝛾_𝑃1≥5), and constant elsewhere. Figure 8(e) is plotted for various values of the reaction/diffusion parameter 𝛾_𝑃1 as abscissa and normalized steady-state current 𝜓 as ordinate. It is observed that the current increases as the value of 𝛾_𝑃1 increases. Also current decreases for all values of the reaction/diffusion parameter 𝛾_𝑃1 and reaches zero when 𝐴≥10.

(a)

(b)

(c)

(d)

(e)

Figure 8 

Dimensionless current 𝜓 versus dimensionless parameters 𝛾𝑃1, (𝛼+𝛽), 𝛾𝑆1, 𝛾𝑆2 and 𝐴 for some fixed values of the other parameters using (17).

7. Conclusion

Mathematical model, describing the steady-state response of an amperometric biosensor based on parallel substrate conversion (catalase-peroxidase biochemical reactions), is presented. A nonlinear time independent coupled system of partial differential equations has been solved analytically using a new Homotopy perturbation method. Analytical expressions for substrate and product concentrations as function of single dimensionless parameter 𝐴 or 𝐵 are derived. The optimum value of 𝐴 or 𝐵 is also reported. The analytical expressions provided are much simpler and more efficient to understand and to predict the behavior of the system. Our analytical solution for substrates and products are compared with numerical solution providing consistent results. Concentration of substrates, products, and current depends upon reaction/diffusion parameters 𝛾_𝑆1, 𝛾_𝑆2, 𝛾_𝑆3, 𝛾_𝑃1, and 𝛾_𝑃2. The response of the steady-state current with respect to all parameters is discussed. For all different concentrations, the analytical results matched well with the simulated result. The close matching of simulated result and the analytical data shows that our proposed solution is able to simulate the dynamic performance of sensing system using the parameters. The analytical results provided here are useful in sensor design, optimization, and prediction of electrode response.

Appendices

A. Approximate Analytical Solution of the Nonlinear Equations (5)–(8) Using the Boundary Conditions (9) by New Homotopy Perturbation Method

In this appendix, we indicated how to find the solution of (5)–(8). We first construct the Homotopy as follows: 𝑑1−𝑝²𝑆_1𝑁𝑑𝑦²−𝛾_𝑆1𝑆_1𝑁−𝛾_𝑆2𝑆_1𝑁𝑆_2𝑁(1)𝛼𝑆_1𝑁(1)+𝛽𝑆_2𝑁(1)+𝑝𝛼𝑆_1𝑁𝑑²𝑆_1𝑁𝑑𝑦²+𝛽𝑆_2𝑁𝑑²𝑆_1𝑁𝑑𝑦²−𝛼𝛾_𝑆1𝑆²_1𝑁−𝛽𝛾_𝑆1𝑆_1𝑁𝑆_2𝑁−𝛾_𝑆2𝑆_1𝑁𝑆_2𝑁=0,𝑑1−𝑝²𝑆_2𝑁𝑑𝑦²−𝛾_𝑆3𝑆_1𝑁(1)𝑆_2𝑁𝛼𝑆_1𝑁(1)+𝛽𝑆_2𝑁(1)+𝑝𝛼𝑆_1𝑁𝑑²𝑆_2𝑁𝑑𝑦²+𝛽𝑆_2𝑁𝑑²𝑆_2𝑁𝑑𝑦²−𝛾_𝑆3𝑆_1𝑁𝑆_2𝑁=0,𝑑1−𝑝²𝑃_1𝑁𝑑𝑦²−𝛾_𝑃1𝑆_1𝑁2𝑑+𝑝²𝑃_1𝑁𝑑𝑦²−𝛾_𝑃1𝑆_1𝑁2=0,𝑑1−𝑝²𝑃_2𝑁𝑑𝑦²−𝛾_𝑃2𝑆_1𝑁(1)𝑆_2𝑁𝛼𝑆_1𝑁(1)+𝛽𝑆_2𝑁(1)+𝑝𝛼𝑆_1𝑁𝑑²𝑃_2𝑁𝑑𝑦²+𝛽𝑆_2𝑁𝑑²𝑃_2𝑁𝑑𝑦²−𝛾_𝑃2𝑆_1𝑁𝑆_2𝑁=0. (A.1)(A.2)(A.3)(A.4) The boundary conditions are as follows: when𝑦=0;𝑑𝑆_1𝑁=𝑑𝑦𝑑𝑆_2𝑁=𝑑𝑦𝑑𝑃_2𝑁𝑑𝑦=0;𝑃_1𝑁=0,when𝑦=1;𝑆_1𝑁=𝑆_2𝑁=1,𝑃_1𝑁=𝑃_2𝑁=0. (A.5) The approximate analytical solutions of (A.1) to (A.4) are given by 𝑆_1𝑁=𝑆_1𝑁0+𝑝𝑆_1𝑁1+𝑝²𝑆_1𝑁2𝑆+⋯_2𝑁=𝑆_2𝑁0+𝑝𝑆_2𝑁1+𝑝²𝑆_2𝑁2𝑃+⋯_1𝑁=𝑃_1𝑁0+𝑝𝑃_1𝑁1+𝑝²𝑃_1𝑁2𝑃+⋯_2𝑁=𝑃_2𝑁0+𝑝𝑃_2𝑁1+𝑝²𝑃_2𝑁2+⋯. (A.6) Substituting (A.6) into (A.1)–(A.4) and arranging the coefficients of powers 𝑝, 𝑝⁰:𝑑²𝑆_1𝑁0𝑑𝑦²−𝛾_𝑆1+𝛾_𝑆2𝑆𝛼+𝛽^1𝑁0𝑝=0,₀:𝑑²𝑆_2𝑁0𝑑𝑦²−𝛾_𝑆3𝑆𝛼+𝛽^2𝑁0𝑝=0,₀:𝑑²𝑃_1𝑁𝑑𝑦²+𝛾_𝑃12𝑆_1𝑁0𝑝=0,₀:𝑑²𝑃_2𝑁𝑑𝑦²−𝛾_𝑆3𝑆𝛼+𝛽^2𝑁0=0. (A.7) Solving (A.7) for the boundary conditions (A.5), we can find the following results: 𝑆_1𝑁0𝑦=cosh𝛾_𝑆1𝛼+𝛽+𝛾_𝑆2/𝛼+𝛽𝑦cosh𝛾_𝑆1𝛼+𝛽+𝛾_𝑆2/,𝑆𝛼+𝛽_2𝑁0𝑦=cosh𝛾_𝑆3/𝛼+𝛽𝑦cosh𝛾_𝑆3/,𝑃𝛼+𝛽_1𝑁0𝑦=𝛾_𝑃1𝛼+𝛽−cosh𝛾_𝑆1𝛼+𝛽+𝛾_𝑆2𝑦𝛼+𝛽+𝑦cosh𝛾_𝑆1𝛼+𝛽+𝛾_𝑆2×2𝛾𝛼+𝛽−𝑦+1_𝑆1𝛼+𝛽+𝛾_𝑆2×cosh𝛾_𝑆1𝛼+𝛽+𝛾_𝑆2𝛼+𝛽⁻¹,𝑃_2𝑁0𝑦=−cosh𝛾_𝑆3/𝛼+𝛽𝑦𝛾_𝑃2cosh𝛾_𝑆3/𝛼+𝛽𝛾_𝑆3+𝛾_𝑃2𝛾_𝑆3. (A.8) Since it is a new approach to HPM, we need not to find further iterations. We consider 𝑆_1𝑁0, 𝑆_2𝑁0, 𝑃_1𝑁0, and 𝑃_2𝑁0 as the better approximate solutions. According to the HPM, we can conclude that 𝑆_1𝑁𝑦=lim_𝑝→1𝑆_1𝑁𝑦=𝑆_1𝑁0+𝑆_1𝑁1+⋯≈𝑆_1𝑁0𝑦,𝑆_2𝑁𝑦=lim_𝑝→1𝑆_2𝑁𝑦=𝑆_2𝑁0+𝑆_2𝑁1+⋯≈𝑆_2𝑁0𝑦,𝑃_1𝑁𝑦=lim_𝑝→1𝑃_1𝑁𝑦=𝑃_1𝑁0+𝑃_1𝑁1+𝑃_1𝑁2+⋯≈𝑃_1𝑁0𝑦,𝑃_2𝑁𝑦=lim_𝑝→1𝑃_2𝑁𝑦=𝑃_2𝑁0+𝑃_2𝑁1+𝑃_2𝑁2+⋯≈𝑃_2𝑁0𝑦. (A.9) After putting (A.8) into (A.9), respectively, the final results can be described in (11) to (14) in the text.

B. Matlab Program to Find the Numerical Solution of (5)–(8)

See Algorithm 1.

function pdex9 m = 0; x = linspace(0,1); t  =  linspace(0,10000000000); 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,1) ′ ) %–––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––– 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,3) ′ ) %–––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––– figure plot(x,u4(end,:)) title( ′ u4(x,t) ′ ) xlabel( ′ Distance x ′ ) ylabel( ′ u4(x,4) ′ ) %–––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––– function [c,f,s] = pdex4pde(x,t,u,DuDx) c = [0; 0; 0; 0]; f = [1; 1; 1; 1].* DuDx; r1  =  0.1; r2  =  0.1; r3  =  0.1; r4  =  0.1; r5  =  0.1; a  =  0.1; b  =  0.1; F1  = - r1  *  u(1) - (r2  *  u(1) *  u(2))/(a  *  u(1)+u(2) *  b); F2  = - (r3  *  u(1) *  u(2))/(a  *  u(1)+u(2) *  b); F3  = ((r4/2) *  u(1)); F4  = (r5  *  u(1) *  u(2))/(a  *  u(1)+u(2) *  b); s  = [F1; F2; F3; F4]; % ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––– function u0 = pdex4ic(x)                %create a initial conditions u0 = [1; 1; 0; 0]; % ––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––– function [pl,ql,pr,qr]= pdex4bc(xl,ul,xr,ur,t)      %create a boundary conditions pl = [0; 0; ul(3); 0]; ql = [1; 1; 0; 1]; pr = [ur(1)-1; ur(2)-1; ur(3)-0; ur(4)-0]; qr = [0; 0; 0; 0];

Algorithm 1 

Matlab program to find the numerical solution of (5)–(8).

Nomenclature

𝑆1,𝑒,𝑆2,𝑒:	Substrate concentrations (molcm−3)
𝑒1,𝑒2:	Enzyme concentrations (molcm−3)
𝑃1,𝑒,𝑃2,𝑒:	Product concentrations (molcm−3)
𝐷𝑆1,𝑒,𝐷𝑆2,𝑒,𝐷𝑃1,𝑒,𝐷𝑃2,𝑒:	Diffusion coefficients (cm2sec−1)
𝑘1,𝑘21,𝑘22:	Reaction rate constant (molcm−3)
𝑑𝑒,𝑑𝑑:	Thickness of the enzyme and diffusion layer (cm)
𝑆10,𝑆20:	Bulk concentrations (molcm−3)
𝑥:	Space (cm)
𝑖:	Current (Ampere).

Dimensionless Parameters (Greek Symbols)

𝑆1𝑁,𝑆2𝑁:	Dimensionless substrate concentrations
𝑃1𝑁,𝑃2𝑁:	Dimensionless product concentrations
𝛾𝑆1,𝛾𝑆2,𝛾𝑆3,𝛾𝑃1,𝛾𝑃2:	Dimensionless reaction/diffusion parameters
𝛼,𝛽:	Dimensionless parameters
𝑦:	Dimensionless space
𝐴=(𝛾𝑆1(𝛼+𝛽)+𝛾𝑆2)/(𝛼+𝛽):	Dimensionless parameter
𝐵=𝛾𝑆3/(𝛼+𝛽):	Dimensionless parameter.

Conflict of Interests

The authors declare that they have no conflict of interests regarding the publication of this paper.

Acknowledgment

This work was supported by the CSIR (no. 01 (2442)/10/EMR-II), Government of India.