Mathematical Model for Bio-directional Diffusion of Reactants and Products in the Enzymatic Reaction of Glucose in a Spherical Matrix

In this paper the theoretical model of glucose–oxidaise loaded in chitosan-aliginate microsphere and hydrogen peroxide production is discussed. The glucose and oxygen in the medium diffuse into the microsphere and react, as a catalyst by glucose oxidase, to produce gluconic acid and hydrogen peroxide. The model involves the system of nonlinear nonsteadystate reaction-diffusion equations. Analytical expressions for the concentrations of glucose, oxygen, gluconic acid and hydrogen peroxide are derived from these equations using homotopy perturbation and the reduction of order method. A comparison of the analytical approximation and numerical simulation is also presented. An agreement between analytical expressions and numerical results is observed. The effect of various parameters (glucose concentration in the external solution, particle size, enzyme loading and Michaelis constant etc.) on the concentration of gluconic acid and hydrogen peroxide release is discussed. Sensitivity analysis of parameters is also discussed.


Introduction
Reactive oxygen species (ROS) are chemically reactive molecules containing oxygen. Examples include peroxides, superoxide, hydroxyl radical and singlet oxygen. Reactive oxygen species and cellular oxidants tress have long been associated with cancer. [1] suggest that cancer cells normally produce more ROS than do normal cells. A mathematical model including simultaneous diffusion and enzymatic reaction was developed by [2] Also [3] discuss the theoretical and experimental model of glucose sensitive membrane. [2] have developed a mathematical model to describe a dynamic process of diffusion of reactants, coupled with an enzymatic reaction inside a glucose composite membrane containing an ionic nano particles, glucose oxidase and catalase embedded in a hydrophobic polymer. [4] analyzed the theoretical model describing the process of reaction and diffusion in glucose-responsive composite membranes. Recently, insulin-delivery system contain a glucose membrane has been studied [4][5][6] More recently [7,8] derived the mathematical equation to describe the release kinetics of enzymes, generated hydrogen peroxide from polymeric matrix with spherical geometry. However, to the best of author's knowledge, no general analytical expressions for the concentrations glucose, oxygen, gluconic acid and hydrogen peroxide for all values of reaction/diffusion parameters. The purpose of the present study is to derive simple approximate analytical expression for concentrations profiles of various species inside the GOX-MS for steady and non-steady state conditions. This model can be used to predict the system performance and determine the appropriate combination of material and geometries that can yield useful devices.

Mathematical Formulation of the Problem
Building upon earlier study, [8] presented a concise discussion and derivation of the mass transport nonlinear second-order differential equation which gives the concentration profiles of each species with in the microsphere membrane, is summarized briefly below. The schematic diagram illustrating the structure of a glucose oxidase loaded chitosan-aliginate microsphere and hydrogen oxide production is represented in Fig. 1 [8]. The reaction of glucose oxidation to produce the gluconic acid and hydrogen peroxide in polymeric microspheres can be written as follows: . Schematic diagram illustrating the structure of a glucose oxidase loaded chitosan-aliginate microsphere and hydrogen oxide production (Abdekhodaie, ChengandWu (2015)).
Using mass conservation law by neglecting convection, the nonlinear reaction-diffusion equation in spherical coordinate are given as follows: Where i denotes individual species, e.g. i=g for glucose, i=OX for oxygen, i=a for gluconic acid and i=h for hydrogen peroxide; the stoichiometric coefficients i v , are 1 a v = and 1 h v = [8]; C is the concentration, D is the diffusion coefficient, r is the length parameter, and ℜ is the overall reaction rate which is given by For non steady-state Eq. 2 becomes 2 max 2 2 ( ) Where S is radius of the microsphere, r=0 is the center of the microsphere, and * OX C and * g C are the concentration of oxygen and glucose in the external solution, respectively. We can assume that the diffusion coefficient of glucose, oxygen, gluconic acid and hydrogen peroxide are equal ( ).
Eqs. 4-7 can be written in the following dimensionless form: Here u, v, w and H are the dimensionless concentrations of glucose, oxygen, gluconic acid and hydrogen peroxide respectively, and g and γ γ are the corresponding Thiele modulus. and α β are the dimensionless rate constants. The corresponding initial and boundary conditions(8-10)becomes

Approximate Analytical Expression for the Concentrations Using HPM
Many problems in natural and engineering sciences are modeled by non-linear partial differential equations. In the last two decades with the rapid development of nonlinear science, there has appeared ever increasing interest of physicists and engineers in the analytical techniques for non linear problems. It is well-known, that perturbation methods provide the most versatile tools available in nonlinear analysis of engineering problems [9][10][11][12][13]. The perturbation methods, like the nonlinear analytical techniques, have their own limitations [14,15]. These facts have motivated to suggest alternate techniques, such as variational iteration [16,17] Adomain decomposition [18][19][20][21] and exp function [22]. In order to overcome these drawbacks, combining the standard homotopy and perturbation method, which is called the homotopy perturbation, modifies the homotopy method.
The basic principle of this method is described in Appendix A. We have obtained the analytical expressions of the concentrations of glucose, oxygen and gluconic acid and hydrogen peroxide be solving the non-linear equations (12) to (15) using new approach of homotopy perturbation method [23,24] as follows: Detailed derivations of the dimensionless concentration of glucose, oxygen, gluconic acid and hydrogen peroxide are described in Appendix B-E. During steady state condition (t → ∞ ), from the above results, the analytical expression for concentrations become The concentration of glucose, oxygen, gluconic acid / hydrogen peroxide at the centre of microsphere is given in Table 1.
Gluconic acid or Hydrogen peroxide

Numerical Simulation
The non linear reaction diffusion equations (12) to (15) for the corresponding boundary conditions (16) to (18) are also solved numerically by using Scilab program (AppendixF).
The numerical solutions are compared with our analytical results in Tables 2 to 3 and in Figure 2. The maximum average error between our analytical results and simulation results for the concentration of glucose, oxygen, gluconic acid (or hydrogen peroxide) are1.64%, 0.21% and 1.26 %.

Relation Between Glucose, Oxygen, Gluconic Acid and Hydrogen Peroxide
Equations (20) and (21) represents the relation between glucose and oxygen (Appendix D), glucose and gluconic acid (Appendix E) for all value at time. In the case of steady state, the relation between glucose, oxygen and gluconic acid is

Determination Flux and Release of Kinetics of Hydrogen Peroxide
The flux of hydrogen peroxide from the surface is defined by Fick's first law: The cumulative amount of hydrogen peroxide can be obtained as follows: Cumulative amount of hydrogen peroxide release for unit volume of microsphere When T is large and When k (<0.01) is very small, / 0 t M V → . The eqn (29) represent the new simple analytical expression of cumulative amount of hydrogen peroxide release per unit volume.

Determination of pH Profile Inside the Microspheres
The pH in the presence of gluconic acid is determined by the concentration of buffer ions and gluconic acid in the microsphere.
( ) ( ) Using the equation (21), we obtain the pH in the presence of gluconic acid as follows:

Results and Discussion
Eqns. (19) to (21) represent the new approximate analytical expressions for the concentrations of glucose, oxygen and gluconic acid for all values of the parameters which satisfy the initial and boundary conditions given by the eqns.(10) to (12).
From Fig. 2 (a) and (b), it is observed that the concentrations of glucose and oxygen are depleted at the centre of the microsphere (R=0) as they are consumed by the enzyme reaction. The slope of this decrease in glucose and oxygen increase with the increase in thiele moduli or radius of the microsphere. Since glucose and oxygen together forms gluconic acid at the centre of the microsphere gluconic acid's concentration increases with the increase in the thiele modulus or enzyme concentration.

Influence of Radius of Microsphere over the Concentration of Glucose ( ) u
The radius of chitos analginate microsphere plays a crucial role as the dynamic process involving diffusion of reactants and product is coupled with an enzymatic reaction inside the microsphere. The concentration of glucose depends on the permeability of reactants and products through the membrane. From Fig. 3 it is inferred that when the thiele modulus g or γ γ which depends on the radius of microsphere or enzyme concentration is increased, the concentration of glucose ( ) u decreases. The thicker radius of microsphere is the lower the concentration of glucose. The concentration of glucose drops towards to 0.3 when the value of thiele modulus 70 g γ ≥ or the 100 γ ≥ . The analytical expression of concentration of glucose, oxygen, gluconic acid and hydrogen peroxide versus radius of microsphere is compared with simulation results in Fig. 3, for steady state condition. Satisfactory agreement is noted.

Influence of the Maximal Reaction Velocity (V max ) on the Concentration of Oxygen ( ) v
Increasing the substrate (glucose) concentration in definitely does not increase the rate of an enzyme-catalyzed reaction beyond a certain point. This point is reached when there are enough substrate molecules to completely fill (saturate) the enzyme's active sites. The maximal velocity, or V max , is the rate of the reaction under these conditions. V max reflect show fast the enzyme can catalyze the reaction.
In this enzymatic reaction diffusion process, the maximal reaction velocity max V is proportional to the concentration of the enzyme in the microsphere and the overall kinetics are determined by the maximal reaction rate. From Figs.2 and 3, it is clear that as the reaction diffusion parameter g or γ γ (which depends on the maximal reaction velocity max V or concentration of enzyme) increases, the concentration decreases gradually and becomes zero for higher values of the reaction velocity.

Influence of Concentration of Glucose in the External Solution over the Concentration of Gluconic Acid ( ) w
The concentration of gluconic acid is determined by the concentration of glucose, the permeability of reactants and products through the membrane and the rate of enzymatic reaction.
From Fig.5, it is evident that the concentration of gluconic acid keeps increasing by decreasing the thiele moduli g γ which depends on the initial concentration of glucose * g C .

Influence of Michaelis-Menten Constants on the Concentration of Glucose ( ) u
Michael is constant (K g ) describes the substrate concentration at which half the enzyme's active sites are occupied by substrate. A large K g indicates the need for high substrate concentrations to achieve maximum reaction velocity. A small K g indicates that the enzyme requires only a small amount of substrate to become saturated. Hence, the maximum velocity is reached at relatively low substrate concentrations. Figs. 5-7 shows the plot of the dimensionless concentration of glucose, oxygen and gluconic acid versus dimensionless radius of microsphere for various values of reaction parameter and Michaelis-Menten constant. From these figures it is observed that the concentration glucose and oxygen at the centre of microsphere is decreases when K g or K ox decreases. Increase in Michael is constant value indicates that oxygen also binds with the enzyme.

Influence of Time Parameter on the Glucose, Oxygen, Gluconic Acid
Figs. 4(a)-(b) and Fig. (7) represent the normalized concentration profiles of glucose, oxygen and gluconic acid for various values of time parameter T .From this figures it is inferred that concentration of glucose & oxygen are at centre of the microsphere decreases when time increases. This is due to consumption of glucose and production of hydrogen peroxide at the centre of the microsphere. From Fig. 7(a) it is observed that when time is increases gluconic acid also decreases because of consumption of oxygen. From Fig. 4(b) it is observed that the concentration of glucose, oxygen and gluconic acid for steady state condition at any point of the microsphere is approximately equal to 2 i.e 2 u v w + + ≈ .  The new analytical expression of concentration of glucose and oxygen, at the centre of the microsphere is given in Table-1. Figs. 8 a-b contains plots of glucose and oxygen concentration profiles versus time T. From these graphs, several important observation can be made. The centre of the sphere glucose and oxygen concentration increases from its initial value and reaches the maximum level when 0.1 T = .And then decreases smoothly from the bulk value. A steady state distribution of both glucose and oxygen is achieved when a dimensional time T=0.2. This is due to the constant balance between diffusional supply and enzymatic consumption of oxygen and glucose. In Figs. 9 (a)-(c) all the above discussed results are confirmed.

Influence of Various Parameter on the Release of Hydrogen Peroxide
Figs. 10 shows cumulative release of hydrogen peroxide versus time for various values concentration of glucose in the external solution, different microspheres sizes and varies dimensional parameters k. From this Fig. 10(a) it is observed that the glucose concentration in the external solution increases the hydrogen peroxide release rate. This is due to decrease of internal pH and less hydrogen peroxide permeability. The dimensionless parameter k depends upon thiele modules and Michaelis constants. Thiele modules are directly proportional to maximum reaction rate and radius of the microsphere where as it is inversely proportional to diffusion coefficients and concentration of glucose and oxygen in the external solution. From this Fig. 10(b) it is inferred that as the maximum reaction rate or radius of the microsphere increases the hydrogen peroxide release rate increases. Fig.  10(c), represent cumulative release of hydrogen peroxide for various values radius of microspheres. From the Fig. 10(c), it is observed that the radius of microsphere or particle size decreases, the release rate of hydrogen peroxide increases.
Since when size of the microsphere is small, more glucose and oxygen molecules diffuse. In Fig. 12 the experimental results and analytical results of cumulative release of hydrogen peroxide are compared. Satisfactory agreement is noted. Fig. 11

Differential Sensitivity Analysis of Parameters
One approach to sensitivity analysis is local sensitivity analysis, which is derivative based (numerical or analytical). Mathematically, the sensitivity of the cost function with respect to certain parameters is equal to the partial derivative of the function with respect to those parameters. The term local refers to the fact that all derivatives are taken at a single point. For simple functions, this approach is efficient. We   Figure 13. Sensitivity analysis for evaluating the influence of the parameters in the concentration of (a) glucose u(R,T) using equation (19).(b)oxygen v(R,T) using equation (20).(c) gluconic acid w(R,T)=hydrogen peroxide H(R,T)using equation (21).

Conclusion
A mathematical model of glucose oxidase-loaded microsphere is discussed. This model describes the release kinetics of generated H 2 O 2 from GOX-loaded spherical polymeric matrix.
A system of non-linear reaction diffusion equations for non-steady state conditions depicted in this model has been solved analytically and from the obtained analytical results the time and position dependent concentrations and diffusivity of glucose, oxygen and gluconic acid have been predicted. By comparing the numerical results with analytical results the accuracy has been verified. The influence of various parameters (glucose concentration in the external solution, particle size, enzyme loading, Michaelis constant etc) on the concentration of gluconic acid and hydrogen peroxide release are discussed. Also this theoretical model is useful for optimizing the performance of hydrogen peroxide delivery system and to obtain the parameters required for improving the design of the system.

Appendix A: Basic Concepts of the HPM
The HPM method has overcome the limitations of traditional perturbation methods. It can take full advantage of the traditional perturbation techniques, so a considerable deal of research has been conducted to apply the Homotopy technique to solve various strong non-linear equations [12][13][14][15]. To explain this method, let us consider the following function: With the boundary conditions of 0 ( , ) 0, Where 0 D is a general differential operator, 0 B is a boundary operator, f(r) is a known analytical function and Γ is the boundary of the domain Ω . Generally speaking, the operator 0 D can be divided into a linear part L and nonlinear part N. Eq.(A1) can therefore be written as: By the Homotopy technique, we construct a Homotopy ( , ) : is an embedding parameter, and 0 u is an initial approximation of Eq.(A1) that satisfies the boundary conditions. From Eqs.(A4) and (A5),we have When 0 p = ,Eqs.(A4) and (A5) become linear equations. When 1 p = ,they become non-linear equations. The process of changing p from zero to unity is that of ...
This is the basic idea of the HPM.

Appendix B: Approximate Analytical Solution Eqn (12) using HPM Method
We construct the new homotopy for the eqn (12) as follows (Rajendran and Anitha (2013)): is an embedding parameter. Now assume that the solution of the eqns (12)and (13) In Laplace plane this equation becomes To illustrate the basic concepts of reduction of order, we consider the equation Where P, Q, R are function of R. Using reduction of order, from eq(B1),we have 2 ; Be the general solution of eq(B4). If v is so chosen that Substituting the value of P in the above eq(B11),we get Then eq (B8) reduces to

Appendix C: Inverse Laplace Transform of Eqn (B19) Using Complex Inversion Formula
In this appendix we indicate how equation (B15) may be inverted using the complex inversion formula. If ( ) y s represents the Laplace transform of a function ( ) y τ ,then according to the complex inversion formula we can state that The real number c is chosen such that s c = 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 equation (C1), which is then evaluated using the so-called Bromwichcontour.Thecontourintegralisthenevaluatedusingth eresiduetheoremwhichstatesforanyanalyticfunction ( ) The first residue in equation (C5) (12) and (13), we can obtain the following equation.
By using the boundary conditions in (16) to (18), the boundary condition for N will be Now by applying Laplace transform in the eqn (E3) and using the boundary extension (E4) to (E6) and by using complex inversion formula, the solution of (E3) becomes as ( )