Issue
4open
Volume 2, 2019
Difference & Differential Equations and Applications
Article Number 26
Number of page(s) 12
Section Mathematics - Applied Mathematics
DOI https://doi.org/10.1051/fopen/2019015
Published online 11 July 2019

© Q. Wang & Y. Liu, Published by EDP Sciences, 2019

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Introduction

In order to build a platform for biosensers based on a bi-enzyme electrode a research group at Dublin City University in Ireland have conducted a series of experiments which motivated this reseach work (refer to [1, 2] for details). The biosensor system investigated in this paper consists of two enzymes: glucose oxidase (GOX) and horseradish peroxidase (HRP), they are both immobilised onto an electrode modified with a conducting polymer. The two enzymes have very different kinetic characteristics, the performance of the system will depend on the optimal ratio of the two enzymes. A cascade reaction takes place at the electrode, GOX catalyses the oxidation reaction of glucose to gluconic acid, with production of H2O2. HRP is oxidised by hydrogen peroxide and then subsequently reduced by electrons provided by the electrode, which can be summarised as:

β - D - glucose   + O 2 + H 2 O   GOX     gluconic   acid   + H 2 O 2 , H 2 O 2 +   HRP     Compound   I   + H 2 O , Compound   I   + 2 e - + H +   HRP   + H 2 O . $$ \begin{array}{c}\beta -D-\mathrm{glucose}\enspace +{\mathrm{O}}_2+{\mathrm{H}}_2\mathrm{O}\stackrel{\enspace \mathrm{GOX}\enspace }{\to }\enspace \mathrm{gluconic}\enspace \mathrm{acid}\enspace +{\mathrm{H}}_2{\mathrm{O}}_2,\\ {\mathrm{H}}_2{\mathrm{O}}_2+\enspace \mathrm{HRP}\enspace \to \enspace \mathrm{Compound}\enspace I\enspace +{\mathrm{H}}_2\mathrm{O},\\ \mathrm{Compound}\enspace I\enspace +2{\mathrm{e}}^{-}+{\mathrm{H}}^{+}\to \enspace \mathrm{HRP}\enspace +{\mathrm{H}}_2\mathrm{O}.\end{array} $$

Modelling strategies

A number of simplifying assumptions are made in order to carry out mathematical modelling for these reactions (refer to [3] for details). The standard Michaelis–Menten equation (1) was used to model this biochemical reaction:

E 1 + S 1   k 1 k - 1   C 1 k 2 E 1 + S 2 , E 2 + S 2   k 3 k - 3   C 2 k 4 E 2 + P , $$ \begin{array}{cc}{E}_1+{S}_1\enspace \begin{array}{c}{k}_1\\ \rightleftarrows \\ {k}_{-1}\end{array}\enspace {C}_1\stackrel{{k}_2}{\to }{E}_1+{S}_2,& {E}_2+{S}_2\enspace \begin{array}{c}{k}_3\\ \rightleftarrows \\ {k}_{-3}\end{array}\enspace {C}_2\stackrel{{k}_4}{\to }{E}_2+P,\end{array} $$(1)where E 1 represents glucose oxidase, E 2 represents horseradish peroxidase, S 1 and S 2 represent the two substrates: glucose and hydrogen peroxide, C 1 and C 2 represent the two complexes and P represents the final product. All the k values denote the reactions rate, which are all constant.

Three models of varying complexity for analysing this optimisation problem were discussed. Analytical solutions were presented with a view to expressing the steady-state current as a function of the ratio ζ of the two immobilised enzymes and thus finding its maximum value.

The first model, the “comprehensive model” assumes that the two substrates, S 1 and S 2, are free to diffuse in the solution and consists of two diffusion equations with the relevant nonlinear reaction-type boundary conditions. This model was proposed and solved numerically in [1]. We then jump to the other extreme and ignore all transport phenomena basically reducing the whole problem to the one-point kinetics of the cascade reaction in the second model, the “simplified model”. This leads to a system of ordinary differential equations which is analysed using a combination of dynamical systems methods, perturbation techniques and numerical simulations. This model was also studied previously in Refs. [2, 3]. Finally, the third model, the “intermediate model” is proposed here for the first time and is a compromise between the two situations discussed above, where we allow one substrate (hydrogen peroxide) to diffuse but assume the other substrate (glucose) is only present at the electrode, in comparison with the first model, we obtain a much simpler system.

It is perhaps instructive to give some motivations regarding the choice of these three models and discuss expectations. We expect the diffusion of the first substrate S 1 to the electrode will have no effect on the equilibrium state except by increasing the time it takes to achieve it. So we expect there will be little difference between the first and third model. On the other hand, neglecting the diffusion of the second substrate S 2 (hydrogen peroxide) in the second model is potentially more serious as this assumes that all S 2 generated in the first reaction is immediately available for the second reaction and this could affect the size of the final steady states.

The comprehensive model

This section reviews the model introduced in [1] (where both substrates diffuse). An additional steady-state analysis of the partial differential equations is also presented.

Review of the comprehensive model

In this model, the cascade scheme (1) is modelled by a system of partial differential equations and boundary conditions representing convective and diffusive transport of the two substrates, glucose and hydrogen peroxide, as well as the reaction kinetics of the bi-enzyme electrode. For simplicity, the convective transport is not explicitly modelled and the flow injection is only reflected in the boundary conditions imposed at the top of the diffusion domain, 0 ≤ x ≤ L. We have also assumed that diffusion is one-dimensional and x measures distance from the electrode.

In what follows, we denote the concentrations of all the chemical species mentioned in the cascade scheme (1) by their corresponding lower case letters (e.g., s 1 = [S 1], etc.). The two substrates satisfy the following diffusion equations:

s 1 ( x , t ) t = D 1 2 s 1 ( x , t ) x 2 , 0 x L , t 0 , s 2 ( x , t ) t = D 2 2 s 2 ( x , t ) x 2 , 0 x L ,   t 0 ,   $$ \begin{array}{c}\begin{array}{ccc}\frac{\mathrm{\partial }{s}_1(x,t)}{\mathrm{\partial }t}={D}_1\frac{{\mathrm{\partial }}^2{s}_1(x,t)}{\mathrm{\partial }{x}^2},& 0\le x\le L,& t\ge 0,\end{array}\hspace{0.5em}\\ \begin{array}{ccc}\frac{\mathrm{\partial }{s}_2(x,t)}{\mathrm{\partial }t}={D}_2\frac{{\mathrm{\partial }}^2{s}_2(x,t)}{\mathrm{\partial }{x}^2},& 0\le x\le L,\enspace & t\ge 0,\end{array}\enspace \end{array} $$where D 1 and D 2 denotes the coefficients of diffusion. The top boundary conditions were formed based on the fact that the first substrate S 1 is continuously injected into the solution at a constant speed, and the second substrate S 2 is flushed away constantly, where,

s 1 ( L , t ) = s 0 , t 0 , s 2 ( L , t ) = 0 , t 0 .     $$ \begin{array}{c}\begin{array}{cc}{s}_1(L,t)={s}_0,& t\ge 0,\end{array}\\ \begin{array}{cc}{s}_2(L,t)=0,& t\ge 0.\end{array}\enspace \end{array}\enspace $$The bottom boundary conditions are:

D 1 s 1 ( 0 , t ) x = k 1 e 1 ( t ) s 1 ( 0 , t ) - k - 1 c 1 ( t ) , D 2 s 2 ( 0 , t ) x = k 3 e 2 ( t ) s 2 ( 0 , t ) - k 2 c 1 ( t ) - k - 3 c 2 ( t ) , $$ \begin{array}{c}{D}_1\frac{\mathrm{\partial }{s}_1(0,t)}{\mathrm{\partial }x}={k}_1{e}_1(t){s}_1(0,t)-{k}_{-1}{c}_1(t),\\ {D}_2\frac{\mathrm{\partial }{s}_2(0,t)}{\mathrm{\partial }x}={k}_3{e}_2(t){s}_2(0,t)-{k}_2{c}_1(t)-{k}_{-3}{c}_2(t),\end{array} $$together with,

d e 1 d t = - k 1 e 1 ( t ) s 1 ( 0 , t ) + ( k - 1 + k 2 ) c 1 ( t ) , d e 2 d t = - k 3 e 2 ( t ) s 2 ( 0 , t ) + ( k - 3 + k 4 ) c 2 ( t ) , d c 1 d t = k 1 e 1 ( t ) s 1 ( 0 , t ) - ( k - 1 + k 2 ) c 1 ( t ) , d c 2 d t = k 3 e 2 ( t ) s 2 ( 0 , t ) - ( k - 3 + k 4 ) c 2 ( t ) , d p d t = k 4 c 2 ( t ) . $$ \begin{array}{c}\frac{\mathrm{d}{e}_1}{\mathrm{d}t}=-{k}_1{e}_1(t){s}_1(0,t)+({k}_{-1}+{k}_2){c}_1(t),\\ \frac{\mathrm{d}{e}_2}{\mathrm{d}t}=-{k}_3{e}_2(t){s}_2(0,t)+({k}_{-3}+{k}_4){c}_2(t),\\ \begin{array}{c}\frac{\mathrm{d}{c}_1}{\mathrm{d}t}={k}_1{e}_1(t){s}_1(0,t)-({k}_{-1}+{k}_2){c}_1(t),\\ \frac{\mathrm{d}{c}_2}{\mathrm{d}t}={k}_3{e}_2(t){s}_2(0,t)-({k}_{-3}+{k}_4){c}_2(t),\\ \frac{\mathrm{d}p}{\mathrm{d}t}={k}_4{c}_2(t).\end{array}\end{array} $$

The initial conditions are:

e 1 ( 0 ) = e 1 0 ,   e 2 ( 0 ) = e 2 0 ,   c 1 ( 0 ) = 0 ,   c 2 ( 0 ) = 0 ,   p ( 0 ) = 0 ,   s 2 ( x , 0 ) = 0 , $$ \begin{array}{ccc}{e}_1(0)={e}_1^0,\enspace & {e}_2(0)={e}_2^0,\enspace & \begin{array}{ccc}{c}_1(0)=0,\enspace & {c}_2(0)=0,\enspace & \begin{array}{cc}p(0)=0,\enspace & {s}_2(x,0)=0,\end{array}\end{array}\end{array} $$

s 1 ( x ,   0 ) = { s 0 , if   x = L , 0 , otherwise , $$ {s}_1\left(x,\enspace 0\right)=\left\{\begin{array}{cc}{s}_0,& \mathrm{if}\enspace x=L,\\ 0,& \mathrm{otherwise},\end{array}\right. $$where e 1 0 $ {e}_1^0$, e 2 0 $ {e}_2^0$ and s 0 are constants. We let,

ζ = e 1 0 e 2 0 , $$ \zeta =\frac{{e}_1^0}{{e}_2^0}, $$which implies,

e 1 0 = 1 + ζ ,   e 2 0 = e 1 + ζ , $$ \begin{array}{cc}{e}_1^0=\frac{{e\zeta }}{1+\zeta },\enspace & {e}_2^0=\frac{e}{1+\zeta },\end{array} $$in which e represents the total amount of e 1 and e 2 consumed in the biochemical reaction, it is a constant parameter which can be measured experimentally.

The system is non-dimensionalised by the variables as shown below:

s ¯ 1 ( x ¯ , t ¯ ) = s 1 ( x , t ) s 0 , s ¯ 2 ( x ¯ , t ¯ ) = s 2 ( x , t ) s 0 , e ¯ 1 ( t ¯ ) = e 1 ( t ) e , e ¯ 2 ( t ¯ ) = e 2 ( t ) e ,   c ¯ 1 ( t ¯ ) = c 1 ( t ) e , c ¯ 2 ( t ¯ ) = c 2 ( t ) e ,   p ¯ ( t ¯ ) = p ( t ) e ,   x ¯ = x l ,   t ¯ = t t 0 , $$ \begin{array}{c}\begin{array}{ccc}{\bar{s}}_1(\bar{x},\bar{t})=\frac{{s}_1(x,t)}{{s}_0},& {\bar{s}}_2(\bar{x},\bar{t})=\frac{{s}_2(x,t)}{{s}_0},& \begin{array}{cc}{\bar{e}}_1(\bar{t})=\frac{{e}_1(t)}{e},& {\bar{e}}_2(\bar{t})=\frac{{e}_2(t)}{e},\end{array}\end{array}\enspace \\ \begin{array}{ccc}{\bar{c}}_1(\bar{t})=\frac{{c}_1(t)}{e},& {\bar{c}}_2(\bar{t})=\frac{{c}_2(t)}{e},\enspace & \begin{array}{cc}\bar{p}(\bar{t})=\frac{p(t)}{e},\enspace & \bar{x}=\frac{x}{l},\enspace \bar{t}=\frac{t}{{t}_0},\end{array}\end{array}\end{array}\hspace{0.5em} $$where t 0 = 1/(k 1 s 0). We then obtain the non-dimensional system,

{ s 1 ( x ,   t ) t = D 1 k 1 s 0 l 2 2 s 1 ( x , t ) x 2 ( 2 a ) s 2 ( x ,   t ) t = D 2 k 1 s 0 l 2 2 s 2 ( x , t ) x 2 ( 2 b ) s 1 ( 1 , t ) = 1 ( 2 c ) s 2 ( 1 , t ) = 0 ( 2 d ) s 1 ( 0 ,   t ) x = k 1 el D 1 e 1 ( t ) s 1 ( 0 , t ) - k - 1 el D 1 s 0 c 1 ( t ) ( 2 e ) s 2 ( 0 ,   t ) x = k 3 el D 2 e 2 ( t ) s 2 ( 0 , t ) - k 2 el D 2 s 0 c 1 ( t ) - k - 3 el D 2 s 0 c 2 ( t ) ( 2 f ) d e 1 d t = - e 1 ( t ) s 1 ( 0 , t ) + K m 1 s 0 c 1 ( t ) ( 2 g ) d e 2 d t = - k 3 k 1 e 2 ( t ) s 2 ( 0 , t ) + k - 3 + k 4 k 1 s 0 c 2 ( t ) ( 2 h ) d c 1 d t = e 1 ( t ) s 1 ( 0 , t ) - K m 1 s 0 c 1 ( t ) ( 2 i ) d c 2 d t = k 3 k 1 e 2 ( t ) s 2 ( 0 , t ) - k - 3 + k 4 k 1 s 0 c 2 ( t ) ( 2 j ) d p d t = k 4 k 1 s 0 c 2 ( t ) , ( 2 k ) $$ \left\{\begin{array}{c}\begin{array}{cc}\frac{\mathrm{\partial }{s}_1(x,\enspace t)}{\mathrm{\partial }t}=\frac{{D}_1}{{k}_1{s}_0{l}^2}\frac{{\mathrm{\partial }}^2{s}_1(x,t)}{\mathrm{\partial }{x}^2}& (2\mathrm{a})\end{array}\\ \begin{array}{cc}\frac{\mathrm{\partial }{s}_2(x,\enspace t)}{\mathrm{\partial }t}=\frac{{D}_2}{{k}_1{s}_0{l}^2}\frac{{\mathrm{\partial }}^2{s}_2(x,t)}{\mathrm{\partial }{x}^2}& (2\mathrm{b})\end{array}\\ \begin{array}{c}\begin{array}{cc}{s}_1(1,t)=1& (2\mathrm{c})\end{array}\\ \begin{array}{c}\begin{array}{cc}{s}_2(1,t)=0& (2\mathrm{d})\end{array}\\ \begin{array}{c}\begin{array}{cc}\frac{\mathrm{\partial }{s}_1(0,\enspace t)}{\mathrm{\partial }x}=\frac{{k}_1{el}}{{D}_1}{e}_1(t){s}_1(0,t)-\frac{{k}_{-1}{el}}{{D}_1{s}_0}{c}_1(t)& (2\mathrm{e})\end{array}\\ \begin{array}{c}\begin{array}{cc}\frac{\mathrm{\partial }{s}_2(0,\enspace t)}{\mathrm{\partial }x}=\frac{{k}_3{el}}{{D}_2}{e}_2(t){s}_2(0,t)-\frac{{k}_2{el}}{{D}_2{s}_0}{c}_1(t)-\frac{{k}_{-3}{el}}{{D}_2{s}_0}{c}_2(t)& (2\mathrm{f})\end{array}\\ \begin{array}{c}\begin{array}{cc}\frac{\mathrm{d}{e}_1}{\mathrm{d}t}=-{e}_1(t){s}_1(0,t)+\frac{{K}_m^1}{{s}_0}{c}_1(t)& (2\mathrm{g})\end{array}\\ \begin{array}{c}\begin{array}{cc}\frac{\mathrm{d}{e}_2}{\mathrm{d}t}=-\frac{{k}_3}{{k}_1}{e}_2(t){s}_2(0,t)+\frac{{k}_{-3}+{k}_4}{{k}_1{s}_0}{c}_2(t)& (2\mathrm{h})\end{array}\\ \begin{array}{c}\begin{array}{cc}\frac{\mathrm{d}{c}_1}{\mathrm{d}t}={e}_1(t){s}_1(0,t)-\frac{{K}_m^1}{{s}_0}{c}_1(t)& (2\mathrm{i})\end{array}\\ \begin{array}{c}\begin{array}{cc}\frac{\mathrm{d}{c}_2}{\mathrm{d}t}=\frac{{k}_3}{{k}_1}{e}_2(t){s}_2(0,t)-\frac{{k}_{-3}+{k}_4}{{k}_1{s}_0}{c}_2(t)& (2\mathrm{j})\end{array}\\ \begin{array}{cc}\frac{\mathrm{d}p}{\mathrm{d}t}=\frac{{k}_4}{{k}_1{s}_0}{c}_2(t),& (2\mathrm{k})\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right. $$(2)where the bars were dropped for convenience. An extensive numerical analysis of this system was presented in [1] where the behaviour of the enzyme ratio ζ was studied for different values of the system parameters. The time evolution of k 4 c ¯ 2 ( t ) $ {k}_4{\bar{c}}_2(t)$, which was taken as a measure of the amperometric signal, was calculated and the steady state value, k 4 c 2 * $ {k}_4{c}_2^{\mathrm{*}}$, was recorded as the current value and used for future parameter iterations. Note that the dimensional value of the current is,

I ( t ) d p d t = k 4 c 2 ( t ) = e k 4 c ¯ 2 ( t ) , $$ I(t)\approx \frac{\mathrm{d}p}{\mathrm{d}t}={k}_4{c}_2(t)=e{k}_4{\bar{c}}_2(t), $$(3)hence, k 4 c ¯ 2 ( t ) $ {k}_4{\bar{c}}_2(t)$ can be regarded as a good measure of the amperometric signal, which should not depend on the choice of non-dimensionalisation used in the model.

Steady-state analysis

This part of the paper analyse system (2) as t → ∞. At equilibrium, equation (2a) gives,

s 1 ( x , t ) t = 0 , $$ \frac{\mathrm{\partial }{s}_1(x,t)}{\mathrm{\partial }t}=0, $$which means,

2 s 1 ( x , t ) x 2 = 0 ; $$ \frac{{\mathrm{\partial }}^2{s}_1(x,t)}{\mathrm{\partial }{x}^2}=0; $$(4)then integrating equation (4) twice with respect to x, we get,

s 1 * ( x ) = Ax + B . $$ {s}_1^{\mathrm{*}}(x)={Ax}+B. $$(5)Similarly, from equation (2b), we obtain,

s 2 * ( x ) = Cx + D , $$ {s}_2^{\mathrm{*}}(x)={Cx}+D, $$(6)where s 1 * ( x ) $ {s}_1^{\mathrm{*}}(x)$, s 2 * ( x ) $ {s}_2^{\mathrm{*}}(x)$ denote the equilibrium values of s 1(x, t), s 2(x, t) respectively, and A, B, C, D are constants of integration.

Equation (2c) together with equation (5) gives,

s 1 * ( x ) = ( 1 - B ) x + B ; $$ {s}_1^{\mathrm{*}}(x)=\left(1-B\right)x+B; $$similarly equation (2d) together with equation (6) gives,

s 2 * ( x ) = - Dx + D . $$ {s}_2^{\mathrm{*}}(x)=-{Dx}+D. $$Thus, at x = 0, the system (2) can be reduced to,

{ 1 - B = k 1 el D 1 ( ζ 1 + ζ - c 1 * ) B - k - 1 el D 1 s 0 c 1 * ( 7 a ) - D = k 3 el D 2 ( 1 1 + ζ - c 2 * ) D - k 2 el D 2 s 0 c 1 * - k - 3 el D 2 s 0 c 2 * ( 7 b ) 0 = ( ζ 1 + ζ - c 1 * ) B - K m 1 s 0 c 1 * ( 7 c ) 0 = k 3 k 1 ( 1 1 + ζ - c 2 * ) D - k 4 + k - 3 k 1 s 0 c 2 * . ( 7 d ) $$ \left\{\begin{array}{c}\begin{array}{cc}1-B=\frac{{k}_1{el}}{{D}_1}\left(\frac{\zeta }{1+\zeta }-{c}_1^{\mathrm{*}}\right)B-\frac{{k}_{-1}{el}}{{D}_1{s}_0}{c}_1^{\mathrm{*}}& (7\mathrm{a})\end{array}\\ \begin{array}{cc}-D=\frac{{k}_3{el}}{{D}_2}\left(\frac{1}{1+\zeta }-{c}_2^{\mathrm{*}}\right)D-\frac{{k}_2{el}}{{D}_2{s}_0}{c}_1^{\mathrm{*}}-\frac{{k}_{-3}{el}}{{D}_2{s}_0}{c}_2^{\mathrm{*}}& (7\mathrm{b})\end{array}\\ \genfrac{}{}{0pt}{}{\begin{array}{cc}0=\left(\frac{\zeta }{1+\zeta }-{c}_1^{\mathrm{*}}\right)B-\frac{{K}_m^1}{{s}_0}{c}_1^{\mathrm{*}}& (7\mathrm{c})\end{array}}{\begin{array}{cc}0=\frac{{k}_3}{{k}_1}\left(\frac{1}{1+\zeta }-{c}_2^{\mathrm{*}}\right)D-\frac{{k}_4+{k}_{-3}}{{k}_1{s}_0}{c}_2^{\mathrm{*}}.& (7\mathrm{d})\end{array}}\end{array}\right. $$(7)Here, we have a system of four equations with four unknowns B, D, c 1 * $ {c}_1^{\mathrm{*}}$ and c 2 * $ {c}_2^{\mathrm{*}}$, we then reduce the system to two equations in terms of c 1 * $ {c}_1^{\mathrm{*}}$ and c 2 * $ {c}_2^{\mathrm{*}}$, which denote the equilibrium values of c 1(t), c 2(t) respectively.

{ ( c 1 * ) 2 - ( D 1 s 0 k 2 el + D 1 K m 1 k 2 el + ζ 1 + ζ ) c 1 * + D 1 s 0 k 2 el ζ 1 + ζ = 0 ( 8 a ) ( c 2 * ) 2 - ( k 2 c 1 * k 4 + D 2 K m 2 k 4 el + 1 1 + ζ ) c 2 * + k 2 c 1 * k 4 1 1 + ζ = 0 . ( 8 b ) $$ \left\{\begin{array}{c}\begin{array}{cc}({c}_1^{\mathrm{*}}{)}^2-\left(\frac{{D}_1{s}_0}{{k}_2{el}}+\frac{{D}_1{K}_m^1}{{k}_2{el}}+\frac{\zeta }{1+\zeta }\right){c}_1^{\mathrm{*}}+\frac{{D}_1{s}_0}{{k}_2{el}}\cdot \frac{\zeta }{1+\zeta }=0& (8\mathrm{a})\end{array}\\ \begin{array}{cc}({c}_2^{\mathrm{*}}{)}^2-\left(\frac{{k}_2{c}_1^{\mathrm{*}}}{{k}_4}+\frac{{D}_2{K}_m^2}{{k}_4{el}}+\frac{1}{1+\zeta }\right){c}_2^{\mathrm{*}}+\frac{{k}_2{c}_1^{\mathrm{*}}}{{k}_4}\cdot \frac{1}{1+\zeta }=0.& (8\mathrm{b})\end{array}\end{array}\right. $$(8)System (8) can be easily solved to give explicit formulas for c 1 * $ {c}_1^{\mathrm{*}}$ and c 2 * $ {c}_2^{\mathrm{*}}$ where,

c 1 * = ( D 1 s 0 k 2 el + D 1 K m 1 k 2 el + ζ 1 + ζ ) ± ( D 1 s 0 k 2 el + D 1 K m 1 k 2 el + ζ 1 + ζ ) 2 - 4 D 1 s 0 k 2 el ζ 1 + ζ 2 , $$ {c}_1^{\mathrm{*}}=\frac{\left(\frac{{D}_1{s}_0}{{k}_2{el}}+\frac{{D}_1{K}_m^1}{{k}_2{el}}+\frac{\zeta }{1+\zeta }\right)\pm \sqrt{{\left(\frac{{D}_1{s}_0}{{k}_2{el}}+\frac{{D}_1{K}_m^1}{{k}_2{el}}+\frac{\zeta }{1+\zeta }\right)}^2-\frac{4{D}_1{s}_0}{{k}_2{el}}\cdot \frac{\zeta }{1+\zeta }}}{2}, $$(9)

c 2 * = ( k 2 c 1 * k 4 + D 2 K m 2 k 4 el + 1 1 + ζ ) ± ( k 2 c 1 * k 4 + D 2 K m 2 k 4 el + 1 1 + ζ ) 2 - 4 k 2 c 1 * k 4 1 1 + ζ 2 . $$ {c}_2^{\mathrm{*}}=\frac{\left(\frac{{k}_2{c}_1^{\mathrm{*}}}{{k}_4}+\frac{{D}_2{K}_m^2}{{k}_4{el}}+\frac{1}{1+\zeta }\right)\pm \sqrt{{\left(\frac{{k}_2{c}_1^{\mathrm{*}}}{{k}_4}+\frac{{D}_2{K}_m^2}{{k}_4{el}}+\frac{1}{1+\zeta }\right)}^2-\frac{4{k}_2{c}_1^{\mathrm{*}}}{{k}_4}\cdot \frac{1}{1+\zeta }}}{2}. $$(10)Note that the smaller solution is selected in both these quadratic equations, as we need,

c 1 * e 1 0 = ζ 1 + ζ ; c 2 * e 2 0 = 1 1 + ζ .   $$ \begin{array}{cc}{c}_1^{\mathrm{*}}\le {e}_1^0=\frac{\zeta }{1+\zeta };& {c}_2^{\mathrm{*}}\le {e}_2^0=\frac{1}{1+\zeta }.\end{array}\enspace $$We then plot the current, k 4 c 2 * $ {k}_4{c}_2^{\mathrm{*}}$, as a function of ζ by using MAPLE in Figure 1. We note that, as we chose k 2 = k 4, the optimal ratio ζ* approaches one for large glucose concentrations. Figure 2 shows the current as a function of ζ when the ratio k 4/k 2 is varied.

thumbnail Figure 1

Dependence of current on ζ. s 0 = 1 (red), 5 (blue), 10 (black) and 20 (green) mM, k 1 = 102, k −1 = 10−1, k 2 = 10, k 3 = 102, k −3 = 10−1, k 4 = 10, e 0 = 10−5, l = 2 × 10−4, D 1 = 6.7 × 10−10 and D 2 = 8.8 × 10−10.

thumbnail Figure 2

Dependence of current on ζ. k 4/k 2 = 0.2 (red), 0.5 (blue), 1 (black) and 2 (green), k 1 = 102, k −1 = 10−1, k 2 = 10, k 3 = 102, k −3 = 10−1, k 4 = 10, e 0 = 10−5, l = 2 × 10−4, D 1 = 6.7 × 10−10 and D 2 = 8.8 × 10−10.

Simplified model

This section analyses the simple, one-point model of the cascade reaction which was introduced in [2]. A very detailed analysis of the model is presented in [3], the behaviour of the system for different ζ is discussed in here.

Review of the simple model

The following non-dimentional system of ordinary differential equations were constructed for the one-point model based on the fact that the substrates do not diffuse in the biochemical reaction:

{ d s 1 d t = ε 1 ( k - 1 k 1 s 0 c 1 - ( ζ 1 + ζ - c 1 ) s 1 ) ( 11 a ) d s 2 d t = ε 1 ( k 2 k 1 s 0 c 1 - k 3 k 1 ( 1 1 + ζ - c 2 ) s 2 + k - 3 k 1 s 0 c 2 ) ( 11 b ) d c 1 d t = ( ζ 1 + ζ - c 1 ) s 1 - K m 1 s 0 c 1 ( 11 c ) d c 2 d t = k 3 k 1 ( 1 1 + ζ - c 2 ) s 2 - K s 0 c 2 . ( 11 d ) $$ \left\{\begin{array}{c}\begin{array}{cc}\frac{\mathrm{d}{s}_1}{\mathrm{d}t}={\epsilon }_1\left(\frac{{k}_{-1}}{{k}_1{s}_0}{c}_1-\left(\frac{\zeta }{1+\zeta }-{c}_1\right){s}_1\right)& (11\mathrm{a})\end{array}\\ \begin{array}{cc}\frac{\mathrm{d}{s}_2}{\mathrm{d}t}={\epsilon }_1\left(\frac{{k}_2}{{k}_1{s}_0}{c}_1-\frac{{k}_3}{{k}_1}\left(\frac{1}{1+\zeta }-{c}_2\right){s}_2+\frac{{k}_{-3}}{{k}_1{s}_0}{c}_2\right)& (11\mathrm{b})\end{array}\\ \genfrac{}{}{0pt}{}{\begin{array}{cc}\frac{\mathrm{d}{c}_1}{\mathrm{d}t}=\left(\frac{\zeta }{1+\zeta }-{c}_1\right){s}_1-\frac{{K}_m^1}{{s}_0}{c}_1& (11\mathrm{c})\end{array}}{\begin{array}{cc}\frac{\mathrm{d}{c}_2}{\mathrm{d}t}=\frac{{k}_3}{{k}_1}\left(\frac{1}{1+\zeta }-{c}_2\right){s}_2-\frac{K}{{s}_0}{c}_2.& (11\mathrm{d})\end{array}}\end{array}\right. $$(11)We can further simplify system (11) as follows since the first substrate glucose is supplied constantly at the electrode,

{ d c 1 d t = ζ 1 + ζ - ( 1 + K m 1 s 0 ) c 1 ( 12 a ) d c 2 d t = k 3 k 1 ( 1 1 + ζ - c 2 ) s 2 - K s 0 c 2 ( 12 b ) d s 2 d t = ε 1 ( k 2 k 1 s 0 c 1 - k 3 k 1 ( 1 1 + ζ - c 2 ) s 2 + k - 3 k 1 s 0 c 2 ) , ( 12 c ) $$ \left\{\begin{array}{c}\begin{array}{cc}\frac{\mathrm{d}{c}_1}{\mathrm{d}t}=\frac{\zeta }{1+\zeta }-\left(1+\frac{{K}_m^1}{{s}_0}\right){c}_1& (12\mathrm{a})\end{array}\\ \begin{array}{cc}\frac{\mathrm{d}{c}_2}{\mathrm{d}t}=\frac{{k}_3}{{k}_1}\left(\frac{1}{1+\zeta }-{c}_2\right){s}_2-\frac{K}{{s}_0}{c}_2& (12\mathrm{b})\end{array}\\ \begin{array}{cc}\frac{\mathrm{d}{s}_2}{\mathrm{d}t}={\epsilon }_1\left(\frac{{k}_2}{{k}_1{s}_0}{c}_1-\frac{{k}_3}{{k}_1}\left(\frac{1}{1+\zeta }-{c}_2\right){s}_2+\frac{{k}_{-3}}{{k}_1{s}_0}{c}_2\right),& (12\mathrm{c})\end{array}\end{array}\right. $$with initial conditions c 1 (0) = 0, c 2 (0) = 0 and s 2 (0) = 0, where,

ε 1 = e s 0 ,   K m 1 = k - 1 + k 2 k 1 ,   K = k - 3 + k 4 k 1 . $$ {\epsilon }_1=\frac{e}{{s}_0},\enspace {K}_m^1=\frac{{k}_{-1}+{k}_2}{{k}_1},\enspace K=\frac{{k}_{-3}+{k}_4}{{k}_1}. $$

The global behaviour of system (12) is analysed by using the slow-fast dynamics, slow invariant manifold and the dynamical systems analysis (refer to [3] for details), we concluded system (12) has the following long term behaviours,

lim t c 1 ( t ) = c 1 * , for   all ζ ,   $$ \begin{array}{cc}\underset{t\to \mathrm{\infty }}{\mathrm{lim}}{c}_1(t)={c}_1^{\mathrm{*}},& \mathrm{for}\enspace \mathrm{all}\hspace{1em}\zeta,\end{array}\enspace $$(13)

lim t c 2 ( t ) = { c 2 * = ζ ζ * ( 1 + ζ ) , if   ζ ζ * ( 14 a ) 1 1 + ζ , if   ζ ζ * , ( 14 b ) $$ \underset{t\to \mathrm{\infty }}{\mathrm{lim}}{c}_2(t)=\left\{\begin{array}{c}\begin{array}{cc}\begin{array}{cc}{c}_2^{\mathrm{*}}=\frac{\zeta }{{\zeta }^{\mathrm{*}}(1+\zeta )},& \mathrm{if}\enspace {\zeta }\le {\zeta }^{\mathrm{*}}\end{array}& (14\mathrm{a})\end{array}\\ \begin{array}{cc}\begin{array}{cc}\frac{1}{1+\zeta },& \mathrm{if}\enspace {\zeta }\ge {\zeta }^{\mathrm{*}},\end{array}& (14\mathrm{b})\end{array}\end{array}\right. $$(14)

lim t s 2 ( t ) = { s 2 * , if   ζ < ζ * ( 15 a ) , if   ζ ζ * , ( 15 b ) $$ \underset{t\to \mathrm{\infty }}{\mathrm{lim}}{s}_2(t)=\left\{\begin{array}{c}\begin{array}{ccc}{s}_2^{\mathrm{*}},& \mathrm{if}\enspace {\zeta } < {\zeta }^{\mathrm{*}}& (15\mathrm{a})\end{array}\\ \begin{array}{ccc}\mathrm{\infty },& \mathrm{if}\enspace {\zeta }\ge {\zeta }^{\mathrm{*}},& (15\mathrm{b})\end{array}\end{array}\right. $$(15)where,

c 1 * = ζ ( 1 + ζ ) ( 1 + K m 1 s 0 ) , $$ {c}_1^{\mathrm{*}}=\frac{\zeta }{\left(1+\zeta \right)\left(1+\frac{{K}_m^1}{{s}_0}\right)}, $$(16)

c 2 * = s 2 * ( 1 + ζ ) ( s 2 * + K m 2 s 0 ) , $$ {c}_2^{\mathrm{*}}=\frac{{s}_2^{\mathrm{*}}}{\left(1+\zeta \right)\left({s}_2^{\mathrm{*}}+\frac{{K}_m^2}{{s}_0}\right)}, $$(17)

s 2 * = ζ k 2 K m 2 s 0 ( k 4 ( 1 + K m 1 s 0 ) - ζ k 2 ) , $$ {s}_2^{\mathrm{*}}=\frac{\zeta {k}_2{K}_m^2}{{s}_0\left({k}_4\left(1+\frac{{K}_m^1}{{s}_0}\right)-\zeta {k}_2\right)}, $$(18)

ζ * = k 4 k 2 ( 1 + K m 1 s 0 ) . $$ {\zeta }^{\mathrm{*}}=\frac{{k}_4}{{k}_2}\left(1+\frac{{K}_m^1}{{s}_0}\right). $$(19)From these results, we can see that when ζ < ζ*, the initial amount of glucose oxidase is less than the initial amount of horseradish peroxidase. This leads to the amount of s 2 produced from the first reaction is all consummed somehow in the second reaction and eventually an equilibrium state can be reached; and when ζ ≥ ζ*, the initial amount of glucose oxidase is greater than the initial amount of horseradish peroxidase. The production rate is greater than the comsumption rate of s 2, so the concentration of s 2 will increase indefinitely.

We plot the steady state current, k 4 c 2(∞), as a function of ζ for different values of s 0 (Fig. 3) and k 4/k 2 (Fig. 4). In Figure 3 the curves overlay for ζ values ranging from 1 to 6, and in Figure 4 the curves overlay for ζ values ranging from 0 to 1, also note that from equation (14a) and (14b), that the optimal enzymes ratio is given by,

ζ * = k 4 k 2 ( 1 + K m 1 s 0 ) , $$ {\zeta }^{\mathrm{*}}=\frac{{k}_4}{{k}_2}\left(1+\frac{{K}_m^1}{{s}_0}\right), $$as c 2(∞) achieves its maximum value of 1/(1 + ζ*) when ζ = ζ*. Therefore, we are able to obtain an explicit formula which gives the optimal value of ζ in terms of the system parameters. Note the agreement between the results as shown in Figures 3 and 4 and the model in the previous section.

thumbnail Figure 3

Dependence of current on ζ. s 0 = 0.03 (red), 0.09 (blue), 0.2 (black) and 5 (green) mM, k 1 = 102, k −1 = 10−1, k 2 = 10 and k 4 = 10.

thumbnail Figure 4

Dependence of current on ζ. k 4/k 2 = 0.2 (red), 0.5 (blue), 1 (black) and 2 (green), k 1 = 102, k −1 = 10−1, k 2 = 10 and k 4 = 10.

Intermediate model

In this model, we assume the glucose (s 1) does not diffuse but is present only at the electrode boundary point, and in addition it is constant, i.e., s 1(t) = s 0. (In other words, s 1 is supplied continuously at the reaction site.) The second substrate is free to diffuse throughout the domain at all times during the experiment, which is reflected by the following diffusion equation,

s 2 ( x , t ) t = D 1 2 s 2 ( x , t ) x 2 . $$ \frac{\mathrm{\partial }{s}_2(x,t)}{\mathrm{\partial }t}={D}_1\frac{{\mathrm{\partial }}^2{s}_2(x,t)}{\mathrm{\partial }{x}^2}. $$

At the top layer and the electrode, we have the boundary conditions,

s 2 ( L , t ) = 0 , $$ {s}_2(L,t)=0, $$

D 1 s 2 ( 0 , t ) x = k 3 e 2 ( t ) s 2 ( 0 , t ) - k 2 c 1 ( t ) - k - 3 c 2 ( t ) , $$ {D}_1\frac{\mathrm{\partial }{s}_2(0,t)}{\mathrm{\partial }x}={k}_3{e}_2(t){s}_2(0,t)-{k}_2{c}_1(t)-{k}_{-3}{c}_2(t), $$together with,

d e 1 d t = - k 1 e 1 ( t ) s 1 ( t ) + ( k - 1 + k 2 ) c 1 ( t ) , d e 2 d t = - k 3 e 2 ( t ) s 2 ( 0 , t ) + ( k - 3 + k 4 ) c 2 ( t ) , d c 1 d t = k 1 e 1 ( t ) s 1 ( t ) - ( k - 1 + k 2 ) c 1 ( t ) , d c 2 d t = k 3 e 2 ( t ) s 2 ( 0 , t ) - ( k - 3 + k 4 ) c 2 ( t ) , d p d t = k 4 c 2 ( t ) . $$ \begin{array}{c}\frac{\mathrm{d}{e}_1}{\mathrm{d}t}=-{k}_1{e}_1(t){s}_1(t)+({k}_{-1}+{k}_2){c}_1(t),\\ \begin{array}{c}\frac{\mathrm{d}{e}_2}{\mathrm{d}t}=-{k}_3{e}_2(t){s}_2(0,t)+({k}_{-3}+{k}_4){c}_2(t),\\ \begin{array}{c}\frac{\mathrm{d}{c}_1}{\mathrm{d}t}={k}_1{e}_1(t){s}_1(t)-({k}_{-1}+{k}_2){c}_1(t),\\ \begin{array}{c}\frac{\mathrm{d}{c}_2}{\mathrm{d}t}={k}_3{e}_2(t){s}_2(0,t)-({k}_{-3}+{k}_4){c}_2(t),\\ \frac{\mathrm{d}p}{\mathrm{d}t}={k}_4{c}_2(t).\end{array}\end{array}\end{array}\end{array} $$

The initial conditions are:

e 1 ( 0 ) = e 1 0 , e 2 ( 0 ) = e 2 0 ,   s 1 ( 0 ) = s 0 , s 2 ( x , 0 ) = 0 , c 1 ( 0 ) = 0 , c 2 ( 0 ) = 0 , p ( x , 0 ) = 0 , $$ \begin{array}{c}\begin{array}{ccc}{e}_1(0)={e}_1^0,& {e}_2(0)={e}_2^0,\enspace & \begin{array}{cc}{s}_1(0)={s}_0,& {s}_2(x,0)=0,\end{array}\end{array}\hspace{0.5em}\\ \begin{array}{ccc}{c}_1(0)=0,& {c}_2(0)=0,& p(x,0)=0,\end{array}\end{array} $$and the conservation laws are:

{ e 1 ( t ) + c 1 ( t ) = e 1 0 e 2 ( t ) + c 2 ( t ) = e 2 0 . $$ \left\{\begin{array}{c}{e}_1(t)+{c}_1(t)={e}_1^0\\ {e}_2(t)+{c}_2(t)={e}_2^0.\end{array}\right. $$

The system is non-dimensionalised by using the variables below,

s ¯ 2 ( x ¯ , t ¯ ) = s 2 ( x ,   t ) s 0 ,   e ¯ 1 ( t ¯ ) = e 1 ( t ) e ,   e ¯ 2 ( t ¯ ) = e 2 ( t ) e , c ¯ 1 ( t ¯ ) = c 1 ( t ) e ,   c ¯ 2 ( t ¯ ) = c 2 ( t ) e ,   x ¯ = x l ,   t ¯ = t t 0 , $$ \begin{array}{c}\begin{array}{ccc}{\bar{s}}_2(\bar{x},\bar{t})=\frac{{s}_2(x,\enspace t)}{{s}_0},\enspace & {\bar{e}}_1(\bar{t})=\frac{{e}_1(t)}{e},\enspace & {\bar{e}}_2(\bar{t})=\frac{{e}_2(t)}{e},\end{array}\\ \begin{array}{ccc}{\bar{c}}_1(\bar{t})=\frac{{c}_1(t)}{e},\enspace & {\bar{c}}_2(\bar{t})=\frac{{c}_2(t)}{e},\enspace & \begin{array}{cc}\bar{x}=\frac{x}{l},\enspace & \bar{t}=\frac{t}{{t}_0},\end{array}\end{array}\end{array} $$where t 0 = 1/(k 1 s 0); we then obtain the non-dimensional system,

{ s 2 ( x ,   t ) t = D 1 k 1 s 0 l 2 2 s 2 ( x ,   t ) x 2 ( 20 a ) s 2 ( x , 0 ) = 0 ( 20 b ) s 2 ( 1 , t ) = 0 ( 20 c ) s 2 ( 0 ,   t ) x = η ( s 2 ( 0 , t ) ( 1 1 + ζ - c 2 ( t ) ) - κ c 1 ( t ) - μ c 2 ( t ) ) ( 20 d ) d c 1 ( t ) d t = ζ 1 + ζ - ( 1 + σ 1 ) c 1 ( t ) ( 20 e ) d c 2 ( t ) d t = ρ 2 ( s 2 ( 0 , t ) ( 1 1 + ζ - c 2 ( t ) ) - σ 2 c 2 ( t ) ) , ( 20 f ) $$ \left\{\begin{array}{c}\begin{array}{cc}\frac{\mathrm{\partial }{s}_2(x,\enspace t)}{\mathrm{\partial }t}=\frac{{D}_1}{{k}_1{s}_0{l}^2}\frac{{\mathrm{\partial }}^2{s}_2(x,\enspace t)}{\mathrm{\partial }{x}^2}& (20\mathrm{a})\end{array}\\ \begin{array}{c}\begin{array}{cc}{s}_2(x,0)=0& (20\mathrm{b})\end{array}\\ \begin{array}{cc}{s}_2(1,t)=0& (20\mathrm{c})\end{array}\\ \begin{array}{c}\begin{array}{cc}\frac{\mathrm{\partial }{s}_2(0,\enspace t)}{\mathrm{\partial }x}=\eta \left({s}_2(0,t)\left(\frac{1}{1+\zeta }-{c}_2(t)\right)-\kappa {c}_1(t)-\mu {c}_2(t)\right)& (20\mathrm{d})\end{array}\\ \begin{array}{cc}\frac{\mathrm{d}{c}_1(t)}{\mathrm{d}t}=\frac{\zeta }{1+\zeta }-(1+{\sigma }_1){c}_1(t)& (20\mathrm{e})\end{array}\\ \begin{array}{cc}\frac{\mathrm{d}{c}_2(t)}{\mathrm{d}t}={\rho }_2\left({s}_2(0,t)\left(\frac{1}{1+\zeta }-{c}_2(t)\right)-{\sigma }_2{c}_2(t)\right),& (20\mathrm{f})\end{array}\end{array}\end{array}\end{array}\right. $$(20)with non-dimensional initial conditions,

e 1 ( 0 ) = ζ 1 + ζ ,   e 2 ( 0 ) = 1 1 + ζ , s 2 ( x , 0 ) = 0 ,   c 1 ( 0 ) = 0 ,   c 2 ( 0 ) = 0 ,   $$ \begin{array}{ccc}{e}_1(0)=\frac{\zeta }{1+\zeta },\enspace & {e}_2(0)=\frac{1}{1+\zeta },& \begin{array}{cc}{s}_2(x,0)=0,\enspace & \begin{array}{cc}{c}_1(0)=0,\enspace & {c}_2(0)=0,\end{array}\end{array}\end{array}\enspace $$and conservation laws,

{ e 1 ( t ) + c 1 ( t ) = ζ 1 + ζ e 2 ( t ) + c 2 ( t ) = 1 1 + ζ , $$ \left\{\begin{array}{c}{e}_1(t)+{c}_1(t)=\frac{\zeta }{1+\zeta }\\ {e}_2(t)+{c}_2(t)=\frac{1}{1+\zeta },\end{array}\right. $$where,

ζ = e 1 0 e 2 0 , η = k 3 el D 1 ,   κ = k 2 k 3 s 0 , μ = k - 3 k 3 s 0 , ρ 2 = k 3 k 1 , σ 1 = K m 1 s 0 , σ 2 = K m 2 s 0 . $$ \begin{array}{c}\begin{array}{ccc}\zeta =\frac{{e}_1^0}{{e}_2^0},& \eta =\frac{{k}_3{el}}{{D}_1},& \enspace \kappa =\frac{{k}_2}{{k}_3{s}_0},\end{array}\\ \begin{array}{ccc}\mu =\frac{{k}_{-3}}{{k}_3{s}_0},& {\rho }_2=\frac{{k}_3}{{k}_1},& \begin{array}{cc}{\sigma }_1=\frac{{K}_m^1}{{s}_0},& {\sigma }_2=\frac{{K}_m^2}{{s}_0}.\end{array}\end{array}\end{array} $$

We now carry out a steady-state analysis of system (20). At equilibrium,

s 2 ( x , t ) t = 0 , $$ \frac{\mathrm{\partial }{s}_2(x,t)}{\mathrm{\partial }t}=0, $$which gives,

2 s 2 ( x , t ) x 2 = 0 . $$ \frac{{\mathrm{\partial }}^2{s}_2(x,t)}{\mathrm{\partial }{x}^2}=0. $$(21)

Then by integrating (21) twice, we obtain,

s 2 * ( x ) = Ax + B , $$ {s}_2^{\mathrm{*}}(x)={Ax}+B, $$which gives,

s 2 * ( 0 ) = B , $$ {s}_2^{\mathrm{*}}(0)=B, $$(22)where s 2 * ( x ) $ {s}_2^{\mathrm{*}}(x)$ denotes the equilibrium value of s 2(x, t), A and B are constants of integration. Also, from equation (20c), we obtain A = −B and this condition together with equation (22) yields,

s 2 * ( x ) = B ( 1 - x ) , $$ {s}_2^{\mathrm{*}}(x)=B(1-x), $$which gives,

s 2 * ( 0 ) x = - B . $$ \frac{\mathrm{\partial }{s}_2^{\mathrm{*}}(0)}{\mathrm{\partial }x}=-B. $$(23)From equation (20e), we obtain

c 1 * = ζ ( 1 + ζ ) ( 1 + σ 1 ) , $$ {c}_1^{\mathrm{*}}=\frac{\zeta }{(1+\zeta )(1+{\sigma }_1)}, $$and, from equation (20f), we obtain,

c 2 * = 1 1 + ζ s 2 * ( 0 ) s 2 * ( 0 ) + σ 2 . $$ {c}_2^{\mathrm{*}}=\frac{1}{1+\zeta }\frac{{s}_2^{\mathrm{*}}(0)}{{s}_2^{\mathrm{*}}(0)+{\sigma }_2}. $$Also, from equations (20d) and (23), we get,

- B = η ( s 2 * ( 0 ) ( 1 1 + ζ - c 2 * ) - κ c 1 * - μ c 2 * ) , $$ -B=\eta \left({s}_2^{\mathrm{*}}(0)\left(\frac{1}{1+\zeta }-{c}_2^{\mathrm{*}}\right)-\kappa {c}_1^{\mathrm{*}}-\mu {c}_2^{\mathrm{*}}\right), $$which gives,

( 1 + ζ ) B 2 + ( σ 2 ( 1 + ζ ) + η ( σ 2 - μ ) - κ ζ η 1 + σ 1 ) B - κ ζ η σ 2 1 + σ 1 = 0 , $$ (1+\zeta ){B}^2+\left({\sigma }_2(1+\zeta )+\eta ({\sigma }_2-\mu )-\frac{{\kappa \zeta \eta }}{1+{\sigma }_1}\right)B-\frac{{\kappa \zeta \eta }{\sigma }_2}{1+{\sigma }_1}=0, $$from which, B can be easily obtained as a function of ζ, i.e., B(ζ) (Note that since the above quadratic equation has two real roots of different signs, we choose the positive root). Thus, the equilibrium value for the current is,

I k 4 c 2 * = k 4 1 + ζ B ( ζ ) B ( ζ ) + σ 2 . $$ I\approx {k}_4{c}_2^{\mathrm{*}}=\frac{{k}_4}{1+\zeta }\frac{B(\zeta )}{B(\zeta )+{\sigma }_2}. $$ Figure 5 shows the plot of current ( k 4 c 2 * $ {k}_4{c}_2^{\mathrm{*}}$) on ζ for various values of s 0; if we vary k4/k 2 instead, we obtain the graphs in Figure 6.

thumbnail Figure 5

Dependence of current on ζ. s 0 = 0.03 (red), 0.09 (blue), 0.2 (black) and 5 (green), k 1 = 102, k −1 = 10−1, k 2 = 10, k 3 = 102, k −3 = 10−1, k 4 = 10, e 0 = 10−5, l = 2 × 10−4 and D 1 = 6.7 × 10−10.

thumbnail Figure 6

Dependence of current on ζ. k 4/k 2 = 0.2 (red), 0.5 (blue), 1 (black) and 2 (green), k 1 = 102, k −1 = 10−1, k 2 = 10, k 3 = 102, k −3 = 10−1, k 4 = 10, e 0 = 10−5, l = 2 × 10−4 and D 1 = 6.7 × 10−10.

Summary and comparisons

The work presented here shows the study of a bi-enzyme electrode based on a flow injection analysis, with the aim of finding the optimal ratio of the two enzymes involved. Three different mathematical models were considered: the “comprehensive model” (where we include the diffusion effects of the two substrates: glucose and hydrogen peroxide), the “simplified model” (which concentrated on the kinetics of the two reactions and no transport was taken into account) and the “intermediate model” (which only considered the diffusion of the substrate hydrogen peroxide). As the simplified model consisted of a system of ordinary differential equations, we were able to present a detailed analytical study of its solutions (including an exact formula for the optimal ratio of the two enzymes), unlike in the other two models where the results were mostly numerical.

The dependence of the biosensor response (i.e., the measured amperometric current) as a function of ζ, the ratio of the two enzymes, is again plotted in Figure 7, for different values of the glucose concentration and, in Figure 8, for different values of k 4/k 2. Figures 7a and 8a are plotted based on the experimental data from [1].

thumbnail Figure 7

Dependence of current on ζ for different initial concentrations of s 0. Steady-state analysis of (a) and (b) comprehensive model, (c) simplified model, and (d) intermediate model.

thumbnail Figure 8

Dependence of current on ζ for different values of k 4/k 2. Steady-state analysis of (a) and (b) comprehensive model, (c) simplified model, and (d) intermediate model.

To facilitate the comparison of the three models, we then plot the optimal ζ value (the value which maximises the current) as a function of the glucose concentration. The four resulting curves are as shown in Figure 9.

thumbnail Figure 9

Dependence of optimal enzyme ratio on s 0 (glucose oxidase concentration). (a) Steady-state analysis of intermediate model, (b) numerical analysis of comprehensive model, (c) steady-state analysis of comprehensive model, and (d) steady-state analysis of simplified model.

We note that the simple model and the intermediate model give identical results for the optimal enzyme ratio for all values of initial glucose concentration. The values predicted by the comprehensive model are quite different at low glucose concentrations but, again, identical at high glucose concentrations. Note that, at high glucose concentration, the optimal enzyme ratio approaches the same value regardless of the model used (This value is ζ* = 1 in our graph, as a consequence of choosing k 2 = k 4.). The fact that a given optimal enzyme ratio is achieved at a higher value of glucose concentration in the comprehensive model is quite obvious, since in that model the glucose is diffusing from a distant place (unlike in the other two models where s 0 represents the glucose concentration at the reaction site). Also, at high glucose concentrations, we expect both enzymes to be saturated with the corresponding substrates (i.e., working at maximum capacity) and so increasing the amount of glucose will not make any difference to the biosensor performance. The three models will give the same result in this regime as the transport effects only affect the availability of substrates at the reaction site.

Next, we studied the dependence of the optimal enzyme ratio on a different parameter associated with our chemical system, namely k 4/k 2 which represents the ratio of the catalytic turnover numbers for the two consecutive reactions. The optimal ratio ζ was plotted against k 4/k 2 for the three models and the resulting graphs are shown in Figure 10. Note that, the simple model predicts a linear relationship between ζ* and k 4/k 2, as illustrated by equation (19). The three models seem to give identical results at low values of k 4/k 2, but diverge when the second reaction becomes much faster than the first.

thumbnail Figure 10

Dependence of optimal enzyme ratio on k 4/k 2 ratio. (a) Steady-state analysis of intermediate model, (b) numerical analysis of comprehensive model, (c) steady-state analysis of comprehensive model, and (d) steady-state analysis of simplified model.

In conclusion, parameter regimes which are characterised by high glucose concentrations and low values of k 4/k 2 seem relatively indifferent to the modelling strategy used and so we would recommend the simple model, which is the easiest to analyse. For all the other parameter regions, further understanding of the behaviour of the system is required and our results so far seem to imply that we need to investigate the relationship between diffusion rates and the speeds of the two reactions. An asymptotic analysis based on these parameters will form the subject of a future study.

References

  1. Mackey D, Killard AJ, Ambrosi A, Smyth MR (2007), Optimizing the ratio of horseradish peroxidase and glucose oxidase on a bienzyme electrode: Comparison of a theoretical and experimental approach. Sens Actuators B122, 395–402. [CrossRef] [Google Scholar]
  2. Mackey D, Killard AJ (2008), Optimising design parameters of enzyme-channelling biosensors, in: Mattheij R, et al. (Eds.), Progress in Industrial Mathematics at ECMI 2006, Springer, Berlin, Heidelberg, 853–857. [CrossRef] [Google Scholar]
  3. Wang Q, Liu Y (2017), Mathematical model for optimising bi-enzyme biosensors, in: Differential and Difference Equations with Applications, Springer, Cham, 595–615. [Google Scholar]

Cite this article as: Wang Q & Liu Y 2019. Mathematical models for optimising bi-enzyme biosensors. 4open, 2, 26.

All Figures

thumbnail Figure 1

Dependence of current on ζ. s 0 = 1 (red), 5 (blue), 10 (black) and 20 (green) mM, k 1 = 102, k −1 = 10−1, k 2 = 10, k 3 = 102, k −3 = 10−1, k 4 = 10, e 0 = 10−5, l = 2 × 10−4, D 1 = 6.7 × 10−10 and D 2 = 8.8 × 10−10.

In the text
thumbnail Figure 2

Dependence of current on ζ. k 4/k 2 = 0.2 (red), 0.5 (blue), 1 (black) and 2 (green), k 1 = 102, k −1 = 10−1, k 2 = 10, k 3 = 102, k −3 = 10−1, k 4 = 10, e 0 = 10−5, l = 2 × 10−4, D 1 = 6.7 × 10−10 and D 2 = 8.8 × 10−10.

In the text
thumbnail Figure 3

Dependence of current on ζ. s 0 = 0.03 (red), 0.09 (blue), 0.2 (black) and 5 (green) mM, k 1 = 102, k −1 = 10−1, k 2 = 10 and k 4 = 10.

In the text
thumbnail Figure 4

Dependence of current on ζ. k 4/k 2 = 0.2 (red), 0.5 (blue), 1 (black) and 2 (green), k 1 = 102, k −1 = 10−1, k 2 = 10 and k 4 = 10.

In the text
thumbnail Figure 5

Dependence of current on ζ. s 0 = 0.03 (red), 0.09 (blue), 0.2 (black) and 5 (green), k 1 = 102, k −1 = 10−1, k 2 = 10, k 3 = 102, k −3 = 10−1, k 4 = 10, e 0 = 10−5, l = 2 × 10−4 and D 1 = 6.7 × 10−10.

In the text
thumbnail Figure 6

Dependence of current on ζ. k 4/k 2 = 0.2 (red), 0.5 (blue), 1 (black) and 2 (green), k 1 = 102, k −1 = 10−1, k 2 = 10, k 3 = 102, k −3 = 10−1, k 4 = 10, e 0 = 10−5, l = 2 × 10−4 and D 1 = 6.7 × 10−10.

In the text
thumbnail Figure 7

Dependence of current on ζ for different initial concentrations of s 0. Steady-state analysis of (a) and (b) comprehensive model, (c) simplified model, and (d) intermediate model.

In the text
thumbnail Figure 8

Dependence of current on ζ for different values of k 4/k 2. Steady-state analysis of (a) and (b) comprehensive model, (c) simplified model, and (d) intermediate model.

In the text
thumbnail Figure 9

Dependence of optimal enzyme ratio on s 0 (glucose oxidase concentration). (a) Steady-state analysis of intermediate model, (b) numerical analysis of comprehensive model, (c) steady-state analysis of comprehensive model, and (d) steady-state analysis of simplified model.

In the text
thumbnail Figure 10

Dependence of optimal enzyme ratio on k 4/k 2 ratio. (a) Steady-state analysis of intermediate model, (b) numerical analysis of comprehensive model, (c) steady-state analysis of comprehensive model, and (d) steady-state analysis of simplified model.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.