EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 3, Article Number 6203 ISSN 1307-5543 – ejpam.com Published by New York Business Global Mathematical Modeling and Analysis of Immobilized Glucose Dehydrogenase Glucose and Oxygen-Driven Reactions in Spherical Micro-reactors Daniel S1,2, Mallikarjuna M1,3, Senthamarai R1,∗ 1 Department of Mathematics, College of Engineering and Technology, SRM Institute of Science and Technology, Kattankulathur – 603 203, Tamilnadu, India 2 PG & Research Department of Mathematics, Cardamom Planters’ Association College, Bodinayakanur-625 582, Tamil Nadu, India 3 Department of Mathematics, Karpaga Vinayaga College of Engineering and Technology, Madhuranthagam, Tamil Nadu 603308, India Abstract. In this article, a batch stirred tank reactor comprising an array of uniform spherical porous micro-reactor (MR) immobilized with nonspecific glucose dehydrogenase and an oxygen- reducing enzyme is modeled mathematically. This model is a steady-state system of non-linear reaction diffusion terms related to the non-Michaelis-Menten kinetics of an enzymatic reaction. We provide approximate analytical expression of the non-linear reaction diffusion equations of substrate, product and oxygen concentrations by utilizing the Adomian Decomposition Method (ADM) and Taylor Series Method (TSM). The obtained semi-analytical results are compared with the numerical simulations, satisfactory results are obtained. Semi-analytical expression for the effectiveness factor of the system is also provided. Using this expression, the factors influencing the effective operation of the system are discussed. 2020 Mathematics Subject Classifications: 34A34, 34A45, 34B15, 34B60, 92C45 Key Words and Phrases: Micro-reactor, Analytical solution, Taylor series Method, Adomian Decomposition Method, Adomian polynomials 1. Introduction Micro-reactors (MR) with immobilized enzymes enable selective substrate conversions, efficient utilization of small sample and reagent volumes, cost reduction, shorter process times, and a compact system design [1, 2]. Enzyme immobilization may lead to some loss of activity; however, it often enhances stability by confining the enzyme within support materials [3, 4]. For the process of integration of a micro-reactor into a stirred tank reactor, ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i3.6203 Email addresses: daniel.s@cpacollege.org (Daniel S), mallikarjuna363@gmail.com (M. Mallikarjuna) senthamr@srmist.edu.in (Senthamarai R) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 2 of 26 the support material must withstand shear forces [5]. Porous silica-based materials serve as excellent supports due to their tunable properties, enabling high enzyme loading [1, 6, 7]. The development and enhancement of highly efficient and productive biological pro- cesses require the assessment and consideration of various physical and biochemical char- acteristics [1, 8]. In numerous bioreactors, both liquid and solid phases coexist, making mass transfer a crucial factor to consider [8, 9]. Pore and particle size play a crucial role in determining the total surface area, directly influencing enzyme binding capacity [10]. One approach to optimizing micro-reactor characteristics is through extensive physical experi- mentation. Alternatively, advanced computational modeling techniques can be employed to simulate and analyze bioreactor processes for further improvement [11]. Recently, various mathematical models representing the batch reactors have been cre- ated. Sánchez et al. [12] compared the continuous-flow activated sludge and sequencing batch reactor processes using math to see how well they remove carbon and nitrogen while looking at different operating conditions. Rivera et al. [13] studied the use of computa- tional fluid dynamics (CFD) for modeling electrochemical reactors to enhance the existing technologies and develop new applications. It highlights CDF-driven optimizations in reactor design, electrode configuration, fluid distribution, and operational conditions for various industrial and environmental processes. Harker et al. [14] presented a determinis- tic method for calculating key parameters (k2, k3) in anaerobic digestion (AD) using the AM2 model, with a focus on volatile fatty acid concentration dynamics. Their approach demonstrated high accuracy, with global error below 2.37%, ensuring reliable parameter estimation. Obtaining an analytical solution for a non-linear differential equation is the constant problem faced by the researchers. Various iterative methods have been utilized recently for obtaining the semi-analytical solutions for the non-linear differential equations, which can be utilized for better understanding of the physical system. Mallikarjuna and Sen- thamarai [15] utilized the TSM and ADM for obtaining the semi-analytical solutions of non-linear differential equations of substrate and product concentrations arising in amper- ometric biosensors in steady-state conditions, and for non-steady-state conditions, they utilized the Laplace homotopy perturbation method [16]. Recently, they also utilized the ADM for analyzing the steady-state non-linear model of the batch reactor performance for the synthesis of cephalexin [17]. Sivakumar et al. [18, 19] provided the semi-analytical expressions for the non-steady-state immobilized enzyme systems with and without the external mass transfer resistance by utilizing the Laplace homotopy perturbation method. WM Al-Hayani [20] presented the combination of the Laplace transform and the homotopy perturbation method for solving the sine-Gordon equation. The variational iteration method was utilized by Senthamarai and Saibhavani [21] for providing the semi-analytical solution for the non-linear differential equation arising in the immobilized enzyme reactor. WM Al-Hayani and MT Younis [22] provided an ap- S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 3 of 26 proximate analytical method for Green’s function by utilizing the homotopy perturbation method. They also provided an approximate analytical solution for the linear and non- linear parabolic and hyperbolic partial differential equations using the same methodology [23]. El-Ajou et al. [24] introduced the Laplace Residual Power series method to obtain analytical solutions for various ordinary differential equations, including singular, non- singular, linear, and non-linear cases. Olayiwola et al. [25] developed a mathematical model to evaluate the impact of antiviral drugs and monoclonal antibodies on COVID-19 transmission, analyzing stability, sensitivity, and solution uniqueness using the Laplace- ADM and real-world data. Awonusika [26] applied the ADM to obtain analytical solutions for a class of Lane-Emden equations involving Jacobi functions, presenting power series and closed-form solutions for specific cases. Comparisons with existing results confirm the method’s accuracy and reliability in solving these equations across various symmetric and hyperbolic spaces. Mahmood et al. [27] applied the ADM for weakly non-linear systems of harmonic oscillators with separate resonance passages. Waleed Al-Hayani [28] utilized the combination of ADM with the Green’s function for solving the linear and non-linear boundary value problem of twelfth order. Usman et al. [29] utilized the Natural transform and ADM for solving the 2-dimensional heat equation which is a non-homogeneous PDE. Silambuselvi et al. [30] studied the mathematical model of facilitated diffusion process in liquid membrane by utilizing the ADM. They provided semi analytical results for the so- lute, ligand and solute ligand complex and for the facilitator factor of the system. Mebarek et al. [31] utilized the ADM for solving the non-linear differential equations representing the boron diffusion during the boronizing process. To the best of the authors’ knowledge, no semi-analytical solution currently exists for the mathematical model of a microreactor incorporating glucose dehydrogenase and lac- case enzymes immobilized within carbon nanotubes under steady-state conditions. This study aims to bridge that gap by providing a comprehensive analysis of the system, offer- ing valuable insights into the optimization of enzyme-based micro-reactors. In this research article, a mathematical model of a micro-reactor with the bienzyme system of glucose dehydrogenase and lactase immobilized into the carbon nanotubes is analyzed. The major contribution of this research is that the non-linear term in the reaction-diffusion equation is based on the non-Michaelis-Menten type enzymatic kinetics. The semi-analytical expression for the non-linear reaction diffusion equation is obtained by utilizing the TSM and ADM. Obtained analytical results are compared with the nu- merical simulation, and it is seen that the results correlate satisfactorily for all the values of parameters appearing in the system. The analytical expression for the effective work- ing of the micro-reactor is also obtained and vanished. The effect of parameters on the concentrations and effective working of the system are analyzed by utilizing the obtained analytical results. This research helps to understand the theoretical and practical optimization of the bioenzyme systems immobilized on the carbon nanotubes. The semi-analytical solutions S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 4 of 26 obtained by utilizing the ADM and TSM offer a valuable tool for predicting substrate and product concentrations under varying parametric conditions, which is helpful in optimiz- ing the micro-reactor. These findings have a broad range of applications in biomedical engineering for developing more enhanced biosensors and bioreactors, pharmaceuticals for drug synthesis using enzymes, and wastewater treatment. This study provides a strong analytical foundation for improvising enzyme-based reaction systems in diverse biotech- nological applications. This research work is articulated as in Section 2, the mathematical model is given; in Section 3, the semi-analytical solution of the non-linear reaction diffusion equation is depicted; Section 4 depicts the validation of the analytical results; Section 5 depicts the major results of the article; and the conclusion is given in Section 6. 2. Mathematical Formulation of the Problem The subsequent biochemical reactions that occur in the micro-reactor are taken into ac- count for the modelling. Assuming all the following equations satisfy the Michaelis–Menten kinetics. L+ E1OX k1−→ P + E1red, (1) E1red k2−→ E1OX + 2e−, (2) E2OX + 4e− k3−→ E2red, (3) E2red +O2 k4−→ E2OX + 2H2O, (4) where L, O2, and H2O denote lactose, oxygen, and water. The oxidized and reduced enzymes for GDH are E1OX and E1red, E2red stands for lactose, and P is the reaction product. Also, k1 and k4 stand for rate reaction constants of enzyme reactions, and k2 and k3 stand for rate reaction constants of enzyme transfer. Outside of the MR, no enzymatic or energy transfer reactions occur. The Nernst diffusion layer, a thin spherical shell adjacent to the MR surface, remains stagnant over time, assuming that the solution is constantly stirred and the Nernst approach is applied. Therefore, only the mass transport by diffusion occurs in the shell that encircles the MR. When the solution is outside the diffusion shell, it is moving, and the amounts of all the soluble substances are the same all around outside the shell. The schematic representation of immobilized glucose spherical micro-reactor is represented in Fig. 1. S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 5 of 26 Figure 1: Schematic representation of immobilized glucose spherical micro-reactor Nomenclature Symbol Meaning Lm Concentration of Lactose Pm Concentration of Product Om Concentration of Oxygen Lb, Pb, Ob Bulk concentrations r Spherical Co-ordinate R Thickness of Spherical Co-ordinate VL, VO Volumetric Reaction Rate KL, KO Michaelis-Menten Constants A Dimensionless Concentration of Lactose C Dimensionless Concentration of Product B Dimensionless Concentration of Oxygen χ Dimensionless Spherical Co-ordinate ϕ2 1, ϕ2 2, ϕ2 3 Dimensionless Diffusion Parameters α, β, γ, σ Dimensionless Saturation Parameters η Effectiveness Factor ADM Adomian Decomposition Method TSM Taylor Series Method MR Micro-Reactor EF Effectiveness Factor Based on the reaction inside and outside the MR, the steady-state mathematical model is formed. Lm, Pm, Om represents the concentration of lactose, product, and oxygen, respectively. S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 6 of 26 The following are the dimensional steady-state differential equations [11] Dl,m ( d2Lm dr2 + 2 r dLm dr ) = V, (5) Dp,m ( d2Pm dr2 + 2 r dPm dr ) = −V, (6) Do,m ( d2Om dr2 + 2 r dOm dr ) = V 2 , (7) where Pm(r) depicts the product concentration of the MR, Lm(r) and Om(r) are the volumetric concentrations of the lactose and oxygen, respectively, and r is the spherical coordinate of the MR, where V = 2VLVOLmOm KOLmVL + 2KLOmVO + LmOmVL + 2LmOmVO , (8) due to the symmetry zero-flux boundary conditions are defined on the center of the MR DS,m dSm dr ∣∣∣∣ r=0 = 0, S = L, P, O. (9) On the boundary between the diffusion and convective shells (r = R1), the continuity of the concentrations is required Sd(R1) = Sb, S = L, P, O, (10) for the ease of the solution we utilize the following non-dimensional parameters for non- dimensionalizing the governing equations: A(χ) = Lm Lb , B(χ) = Om Ob , C(χ) = Pm Pb , χ = r R1 , ϕ2 1 = R2VLV0 Dl,m , ϕ2 2 = R2VLV0 Dp,m , ϕ2 3 = R2VLV0 Do,m , α = VLKOLb, β = ObKLVO, γ = ObLbVL, σ = VOLbOb, (11) the non-dimensional governing equation for the lactose, product and oxygen concentrations is obtained as follows: d2A(χ) dχ2 + 2 χ dA(χ) dχ = ϕ2 1 2A(χ)B(χ) αA(χ) + 2βB(χ) + γA(χ)B(χ) + 2σA(χ)B(χ) , (12) d2C(χ) dχ2 + 2 χ dC(χ) dχ = −ϕ2 2 2A(χ)B(χ) αA(χ) + 2βB(χ) + γA(χ)B(χ) + 2σA(χ)B(χ) , (13) d2B(χ) dχ2 + 2 χ dB(χ) dχ = ϕ2 3 A(χ)B(χ) αA(χ) + 2βB(χ) + γA(χ)B(χ) + 2σA(χ)B(χ) , (14) S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 7 of 26 with boundary conditions, when χ = 0, dN dχ ∣∣∣∣ χ=0 = 0, N = A, C, B, (15) when χ = 1, N(1) = 1, N = A, C, B, (16) The effectiveness factor of the system is given as η = 3(α+ β + γ + σ) 2ϕ2 2 dC dχ ∣∣∣∣ χ=1 , (17) 3. Semi-Analytical Expressions of Concentration of Lactose, Product and Oxygen Finding an exact solution to a non-linear differential equation is highly challenging. Solving a system of non-linear differential equations remains one of the most complex problems faced by researchers today. In this paper we use two different methods to solve the non-linear equations, the first one was Taylor’s series is an infinite accumulation of terms that are expressed on a single point in terms of the derivatives of a function. Without the use of linearization, pertur- bation, or transformation, the analytical solution is generated by the TSM in the form of highly convergent infinite power series with instantly computable terms. The second method we used is the ADM. It was primarily done by George Adomian in the 1980s [32–34]. It generates results without requiring linearization, perturbation, or discretization of the system. It maintains the accuracy of numerical solutions while significantly reducing computational effort. 3.1. Semi-analytical expressions utilizing Taylor Series Method We use the TSM (see Appendix A) to obtain the semi-analytical solution for the governing Eqs. (12) - (16). S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 8 of 26 Solving the Eqs. (12) - (14) for boundary conditions (15) and (16), we get A(χ) = 1 + l(χ− 1) + (χ− 1)2 2! ( −2l + 2ϕ2 1 α+ 2β + γ + 2σ ) + (χ− 1)3 3!(α+ 2β + γ + 2σ)2 (( 6α2 + (24β + 12γ + 24σ)α+ 24β2 + (4ϕ2 1 + 24γ + 48σ)β +6(γ + 2σ)2 ) l + 2ϕ2 1 ((n− 2)α− 4β − 2γ − 4σ) ) + (χ− 1)4 4!(α+ 2β + γ + 2σ)3  − 8(γ + 2σ + α)βϕ2 1l 2 + ( 16β ((n− 1)α− 2β − γ − 2σ)ϕ2 1 −192 (γ 2 + σ + α 2 + β )3) l − 8ϕ2 1 ( −βϕ2 1 + (n− 2)α2+( (n2 + 2n− 8)β + ( n2 2 + n− 4 ) γ + (n2 + 2n− 8)σ − ϕ2 3 4 ) α −8 (γ 2 + σ + β )2)  + (χ− 1)5 5!(α+ 2β + γ + 2σ)4  16β ( (−3l + 7n 2 − 2)α+ (l − 4)β − 3 ( l + 2 3 ) (γ + 2σ) ) ϕ4 1 + ( (40n− 80)α3 + ( (24l3 + (−48n+ 64)l2 + (24n2 − 128n +80)l + 64n2 + 160n− 480)β + (32n2 + 80n− 240)γ +(64n2 + 160n− 480)σ + 2ϕ2 3(n− 4) ) α2 + ( ((48n+ 128)l2 +(−96n2 − 256n+ 320)l + 48n3 + 128n2 + 160n− 960)β2 + ( (48γ + 96σ)l3 − 48(n− 8 3 )(γ + 2σ)l2 + ( (−48n2 − 128n +160)γ + (−96n2 − 256n+ 320)σ + 28ϕ2 3 ) l + (48n3 + 128n2 +160n− 960)γ + (96n3 + 256n2 + 320n− 1920)σ −24ϕ2 3(n+ 2 3 ) ) β + 12 ( (n3 + 8 3 n2 + 10 3 n− 20)γ + ( 2n3 + 16 3 n2 + 20 3 n− 40 ) σ − ϕ2 3(n+ 2 3 ) ) (γ + 2σ) ) α +(320l − 640)β3 + 128 ( l2 + 5 2 l − 15 2 ) (γ + 2σ)β2 +24(l3 + 8 3 l2 + 10 3 l − 20)(γ + 2σ)2β − 80(γ + 2σ)3 ) ϕ2 1 + 1920 (γ 2 + σ + α 2 + β )4 l  (18) S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 9 of 26 C(χ) = 1 +m(χ− 1)− (χ− 1)2 2! ( 2m+ 2ϕ2 2 α+ 2β + γ + 2σ ) + (χ− 1)3 3!(α+ 2β + γ + 2σ)2 ( 24 (γ 2 + σ + α 2 + β )2 m− 4 (n 2 − 1 ) α+ (l − 2)β − γ − 2σ)ϕ2 2 ) + (χ− 1)4 4!(α+ 2β + γ + 2σ)3  ( (8n− 16)α2 + ((8l2 + (−16n+ 16)l + 8n2 + 16n− 64)β +(4n2 + 8n− 32)γ + (8n2 + 16n− 64)σ − 2ϕ2 3)α +(32l − 64)β2 + ((8l2 + 16l − 64)γ + (16l2 + 32l − 128)σ −8ϕ2 1)β − 16(γ + 2σ)2 ) ϕ2 2 − 192 (γ 2 + σ + α 2 + β )3 m  + (χ− 1)5 5!(α+ 2β + γ + 2σ)4  ( (−40n+ 80)α3 + ( ( −24l3 + (48n− 64)l2 + (−24n2 + 128n −80)l − 64n2 − 160n+ 480 ) β + (−32n2 − 80n+ 240)γ +(−64n2 − 160n+ 480)σ − 2ϕ2 3(n− 4))α2 + ( ((−48n− 128)l2 +(96n2 + 256n− 320)l − 48n3 − 128n2 − 160n+ 960)β2 +((−48l3 + (48n− 128)l2 + (48n2 + 128n− 160)l − 48n3 −128n2 − 160n+ 960)γ + (−96l3 + (96n− 256)l2 +(96n2 + 256n− 320)l − 96n3 − 256n2 − 320n+ 1920)σ +(48ϕ2 1 − 28ϕ2 3)l + (−56ϕ2 1 + 24ϕ2 3)n+ 16ϕ2 3 + 32ϕ2 1)β +12 (( −n3 − 8 3 n2 − 10 3 n+ 20 ) γ + (−2n3 − 16 3 n2 − 20 3 n +40)σ + ϕ2 3(n+ 2 3 ) ) (γ + 2σ) ) α+ (−320l + 640)β3 +((−128l2 − 320l + 960)γ + (−256l2 − 640l + 1920)σ −16ϕ2 1(l − 4))β2 − 24 (( l3 + 8 3 l2 + 10 3 l − 20 ) γ + ( 2l3 + 16 3 l2 + 20 3 l − 40 ) σ − 2(l + 2 3 )ϕ2 1 ) (γ + 2σ)β +80(γ + 2σ)3 ) ϕ2 2 + 1920 (γ 2 + σ + α 2 + β )4 m  (19) S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 10 of 26 B(χ) = 1 + n(χ− 1) + (χ− 1)2 2! ( −2n+ ϕ2 3 α+ 2β + γ + 2σ ) + (χ− 1)3 3!(α+ 2β + γ + 2σ)2 ( 6α2 + (ϕ2 3 + 24β + 12γ + 24σ)α+24(β + γ 2 + σ)2)n+ 2ϕ2 3(−α+ (l − 2)β − γ − 2σ) ) + (χ− 1)4 4!(α+ 2β + γ + 2σ)3  − 4αϕ2 3(β + γ 2 + σ)n2 + (8(−α 2 + (l − 1)β − γ 2 − σ)αϕ2 3− 192( γ 2 + σ + α 2 + β)3)n− 4ϕ2 3 ( −αϕ2 3 4 − 2α2 +((l2 + 2l − 8)β − 4γ − 8σ)α+ (4l − 8)β2 + ((l2 + 2l − 8)γ +(2l2 + 4l − 16)σ − ϕ2 1)β − 2(γ + 2σ)2 )  + (χ− 1)5 5!(α+ 2β + γ + 2σ)4  14 (( n 14 − 2 7 ) α+ ( l − 6n 7 − 4 7 ) β − 3(n+ 2 3 )(γ + 2σ) 7 ) αϕ4 3 + ( (20n− 40)α3 + ( ((12l + 32)n2 + (−24l2 − 64l + 80)n+ 12l3 + 32l2 + 40l − 240)β + 16(n2 + 5 2 n− 15 2 )(γ + 2σ) ) α2 + (( 24n3 + (−48l + 64)n2 + (24l2 − 128l + 80)n+ 64l2 + 160l −480)β2 + ( (24γ + 48σ)n3 − 24(γ + 2σ)(l − 8 3 )n2 + ((−24l2 −64l + 80)γ + (−48l2 − 128l + 160)σ + 28ϕ2 1)n+ (24l3 + 64l2 +80l − 480)γ + (48l3 + 128l2 + 160l − 960)σ − 24ϕ2 1 ( l + 2 3 )) β +6(n3 + 8 3 n2 + 10 3 n− 20)(γ + 2σ)2 ) α+ (160l − 320)β3 + ( (64l2 + 160l − 480)γ + (128l2 + 320l − 960)σ +8ϕ2 1(l − 4) ) β2 + 12(γ + 2σ) ( (l3 + 8 3 l2 + 10 3 l − 20)γ +(2l3 + 16 3 l2 + 20 3 l − 40)σ − 2ϕ2 1(l + 2 3 ) ) β − 40(γ + 2σ)3 ) ϕ2 3 + 1920 (γ 2 + σ + α 2 + β )4 n  (20) The constants l,m, n can be determined from Eqs. (18) - (20) for different values of the parameters α, β, γ, σ ϕ2 1, ϕ 2 2 and ϕ2 3. The effectiveness factor is obtained as η = 3 (α+ β + γ + σ) 2ϕ2 2 (1 + l) , (21) 3.2. Semi-analytical expressions utilizing Adomian Decomposition Method We use the ADM (See Appendix B) to obtain the semi-analytical solution for the gov- erning Eqs. (12) - (16). S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 11 of 26 We obtain a solution for A(χ), C(χ) and B(χ) A(χ) = 1 + ϕ2 1(χ 2 − 1) 3(α+ 2β + γ + 2σ) + ϕ2 1(αϕ 2 3 + 4βϕ2 1) 6(α+ 2β + γ + 2σ)3 ( χ4 20 − χ2 6 + 7 60 ) , (22) C(χ) = 1 + ϕ2 2(1− χ2) 3(α+ 2β + γ + 2σ) + ϕ2 2(αϕ 2 3 + 4βϕ2 1) 6(α+ 2β + γ + 2σ)3 ( χ4 20 − χ2 6 − 7 60 ) , (23) B(χ) = 1 + ϕ2 3(χ 2 − 1) 6(α+ 2β + γ + 2σ) + ϕ2 3(αϕ 2 3 + 4βϕ2 1) 12(α+ 2β + γ + 2σ)3 ( χ4 20 − χ2 6 + 7 60 ) , (24) effectiveness factor of the system is obtained as η = 3 (α+ β + γ + σ) 2ϕ2 2 ( − 2ϕ2 2 3α+ 3β + 3γ + 3σ − ϕ2 2 ( αϕ2 3 + 4βϕ2 1 ) 45 (α+ 2β + γ + 2σ)3 ) . (25) 4. Validation of Analytical Results The governing system of non-linear differential Eqs. (12) - (14) have been solved an- alytically by utilizing the semi-analytical methods, the TSM for the Eqs. (18) - (20) and the ADM for the (22) - (24). The numerical simulation of the governing equations has been carried out by utilizing the MATLAB software. The analytical results are compared with the simulated numerical results to obtain satisfactory results. Figs. 2 - 4 represents the comparison between the numerical, TSM, and ADM; it is observed from the graphs that numerical and analytical results correlate satisfactorily. Tables 1 - 3 present the rela- tive error between the numerical, TSM, and ADM. The maximum deviation between the numerical and TSM is 0.063%, while the ADM shows a maximum discrepancy of 0.122%. In contrast to the previous work by Baronas et al. [11], where the numerical method is directly applied to the system without transformation, this study introduces a non- dimensionalized mathematical model. This enhances generalizability and allows system- atic analysis of key parameters such as diffusion and saturation. Additionally, the semi- analytical methods provide explicit closed-form expressions that provide improved inter- pretability of system behaviors under varying parametric conditions. The dual approach of utilizing both ADM and TSM offers an analytical depth and feasibility not explored in the earlier studies. The comparative analysis against numerical results confirms both the accuracy and reliability of the semi-analytical approaches. The results not only validate the methods but also demonstrate their applicability for optimizing reactor design and for understanding the influence of parameters on the system more comprehensively than traditional numerical methods alone. The analytical framework can be extended to simi- lar multi-enzyme systems, offering a general modeling approach for a class of biocatalytic micro-reactors. S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 12 of 26 5. Result and Discussion Figs. 2 - 4 represent the concentration profiles of lactose, product and oxygen in the MR. It is observed that these concentrations depend on the diffusion parameters ϕ2 i , i = 1, 2, 3, saturation parameters α, β, γ, and σ, and the dimensionless distance χ. 5.1. Effect of catalyzed enzyme space The dimensionless distance within the immobilized MR has a major effect on the concentrations and the effective working of the batch reactor. From Figs. 2a - 2d, it can be seen that the increase in dimensionless distance increases the concentration of lactose in the MR and attains its maximum when χ = 1. Figs. 3a - 3d represent the dimensionless concentration of product in the MR. We can see that the dimensionless distance has the inverse effect as it increases, the concentration decreases and attains zero when χ = 1. Figs. 4a - 4d depict the concentration of oxygen; it is observed that as the enzymatic reaction takes place, the oxygen is released from the enzyme-immobilized MR. From these figures it is observed that the oxygen increases with the increase of the dimensionless distance. 5.2. Effect of saturation parameters The relationships between bulk concentrations and the kinetic parameters are repre- sented as α, β, γ, and σ. These parameters characterize the extent to which the catalytic kinetics exhibit saturation (n > 1, where n = α, β, γ, σ) or remain unsaturated (n < 1). Figs. 2, 3 and 4 represent the effect of the saturation parameters on the concentrations of lactose, product, and oxygen. From these figures it is observed that the increase in parameter increases the concentration of lactose and oxygen but has an inverse effect on the concentration of product, as the increase in parameter decreases the concentration of product. We can see that the concentration becomes uniform when α ≥ 100 for lactose, product, and oxygen. The same is observed for the rest of the parameters β, γ, and σ. 5.3. Effect of the diffusion module The diffusion module ϕi, i = 1, 2, 3 quantifies the relation between the maximum reaction rates of lactose and oxygen to the rate of the diffusion through the MR. Table. 1 - 3 represents the effect of the diffusion module on the concentrations of lactose, product, and oxygen. From the Tables. 1 - 3 we can see that the decrease in the values of diffusion parameters ϕ1 and ϕ3 decreases the concentration of lactose and oxygen. But in the case of product concentration, the diffusion module ϕ2 2 has the inverse effect: as the values increase, it decreases the concentration of product. 5.4. Effectiveness factor of the system The effectiveness factor (EF) denotes the ratio of the reaction rate inside the pellet to the reaction rate at its outer surface. Its dimensionless metric ranges from 1.8 to 2, serving S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 13 of 26 as the measure for the catalytic efficiency. The effectiveness factor measures the actual reaction rate within the immobilized pellet relative to the reaction rate in the absence of diffusion limitations. When the entire pellet volume reacts at a high rate, η approaches 2, whereas it nears its minimum value when the reaction within the pellet occurs at a slower rate. The semi-analytical expressions of the effectiveness factor of the MR are obtained using the TSM and ADM. From the Figs. 5a - 5d represents the influence of the diffusion modulus ϕ2 2 and saturation parameters on the effectiveness factor of the system. It is observed from the figures that saturation parameters are increasing functions; as their values increase, the effectiveness factor of the system also increases. The increase in diffusion modulus ϕ2 2 decreases the effectiveness of the MR with carbon nanotubes. 0 0.2 0.4 0.6 0.8 1 0.975 0.98 0.985 0.99 0.995 1 1.005 D im e n s io n le s s C o n c e n tr a ti o n o f L a c to s e = 100 =10 =0.5 (a) 0 0.2 0.4 0.6 0.8 1 Dimensionless Distance 0.975 0.98 0.985 0.99 0.995 1 1.005 D im e n s io n le s s C o n c e n tr a ti o n o f L a c to s e = 100 = 10 = 0.5 (b) 0 0.2 0.4 0.6 0.8 1 Dimensionless Distance 0.975 0.98 0.985 0.99 0.995 1 1.005 D im e n s io n le s s c o n c e n tr a ti o n o f L a c to s e = 0.7 = 10 =100 (c) 0 0.2 0.4 0.6 0.8 1 Dimensionles Distance 0.975 0.98 0.985 0.99 0.995 1 1.005 D im e n s io n le s s c o n c e n tr a ti o n o f L a c to s e = 0.5 = 100 = 10 (d) Figure 2: Dimensionless concentration for Lactose A(χ) versus the dimensionless distance χ for various values of the parameters (2a) α, (2b) β, (2c) γ, (2d) σ where ‘· · · ’ denotes the numerical solution, ‘***’ represents the ADM solution, and ‘⋄⋄⋄’ corresponds to the TSM solution. S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 14 of 26 0 0.2 0.4 0.6 0.8 1 Dimensionless Distance 1 1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 1.1 D im e n s io n le s s C o n c e n tr a ti o n o f P r o d u c t =100 = 10 = 0.5 (a) 0 0.2 0.4 0.6 0.8 1 Dimensionless Distance 1 1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 1.1 D im e n s io n le s s C o n c e n tr a ti o n o f P r o d u c t = 0.6 = 10 = 100 (b) 0 0.2 0.4 0.6 0.8 1 Dimensionless Distance 1 1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 1.1 D im e n s io n le s s c o n c e n tr a ti o n o f P r o d u c t = 0.7 =10 = 100 (c) 0 0.2 0.4 0.6 0.8 1 Dimensionless Distance 1 1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 1.1 D im e n s io n le s s c o n c e n tr a ti o n o f P r o d u c t = 100 = 10 = 0.5 (d) Figure 3: Dimensionless concentration for Product C(χ) versus the dimensionless distance χ for various values of parameters (3a) α, (3b) β, (3c) γ, (3d) σ where ‘· · · ’ denotes the numerical solution, ‘***’ represents the ADM solution, and ‘⋄⋄⋄’ corresponds to the TSM solution. S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 15 of 26 T a b le 1: C om p ar is o n o f n u m er ic a l si m u la ti on fo r th e co n ce n tr at io n of la ct os e A (χ ) w it h an al y ti ca l so lu ti on s ob ta in ed u si n g T S M a n d A D M fo r fi x ed p a ra m et er va lu es α = 0. 5, β = 0. 6, γ = 0. 7, σ = 0. 5, ϕ 2 2 = 0. 6, ϕ 2 3 = 0. 7 an d va ry in g th e ϕ 2 1 . ϕ 2 1 = 0 .1 ϕ 2 1 = 0 .5 ϕ 2 1 = 1 χ N u m er ic a l T S M er ro r % A D M er ro r% N u m er ic a l T S M er ro r % A D M er ro r% N u m er ic a l T S M er ro r% A D M er ro r% 0 .0 0 .9 9 9 0 0 0 .9 9 9 0 0 0 .0 0 0 0 0 0 .9 9 9 0 0 0 .0 0 0 0 0 0 .9 7 5 9 0 0 .9 7 5 7 0 0 .0 2 0 4 9 0 .9 7 5 6 0 0 .0 3 0 7 4 0 .9 0 3 6 0 0 .9 0 4 6 0 0 .1 1 0 6 7 0 .9 0 3 3 0 0 .0 3 3 2 0 0 .2 0 .9 9 9 0 4 0 .9 9 9 0 4 0 .0 0 0 0 7 0 .9 9 9 0 4 0 .0 0 0 0 4 0 .9 7 6 8 7 0 .9 7 6 6 7 0 .0 2 0 1 2 0 .9 7 6 5 9 0 .0 2 8 5 2 0 .9 0 7 4 3 0 .9 0 8 3 7 0 .1 0 4 1 0 0 .9 0 7 3 0 0 .0 1 4 4 3 0 .4 0 .9 9 9 1 6 0 .9 9 9 1 6 0 .0 0 0 2 5 0 .9 9 9 1 6 0 .0 0 0 1 6 0 .9 7 9 7 6 0 .9 7 9 5 8 0 .0 1 8 7 9 0 .9 7 9 5 5 0 .0 2 1 7 5 0 .9 1 8 9 3 0 .9 1 9 7 1 0 .0 8 5 0 7 0 .9 1 9 3 0 0 .0 4 0 3 9 0 .6 0 .9 9 9 3 6 0 .9 9 9 3 5 0 .0 0 0 5 4 0 .9 9 9 3 5 0 .0 0 0 3 4 0 .9 8 4 5 8 0 .9 8 4 4 3 0 .0 1 5 8 7 0 .9 8 4 4 8 0 .0 1 0 1 1 0 .9 3 8 1 5 0 .9 3 8 6 8 0 .0 5 6 6 3 0 .9 3 9 3 4 0 .1 2 7 1 4 0 .8 0 .9 9 9 6 3 0 .9 9 9 6 3 0 .0 0 0 8 7 0 .9 9 9 6 3 0 .0 0 0 5 3 0 .9 9 1 3 3 0 .9 9 1 2 3 0 .0 1 0 2 6 0 .9 9 1 4 0 0 .0 0 6 9 5 0 .9 6 5 1 6 0 .9 6 5 3 9 0 .0 2 4 1 1 0 .9 6 7 4 7 0 .2 3 9 9 2 1 .0 0 .9 9 9 9 9 0 .9 9 9 9 8 0 .0 0 1 1 8 0 .9 9 9 9 8 0 .0 0 0 6 7 1 .0 0 0 0 0 1 .0 0 0 0 0 0 .0 0 0 3 8 1 .0 0 0 3 0 0 .0 3 0 2 8 1 .0 0 0 0 5 1 .0 0 0 0 0 0 .0 0 5 1 9 1 .0 0 3 7 7 0 .3 7 2 3 4 0 .0 0 0 4 9 0 .0 0 0 2 9 0 .0 1 4 3 2 0 .0 0 8 9 8 0 .0 6 2 5 7 0 .1 2 2 0 3 S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 16 of 26 T a b le 2: C o m p a ri so n o f n u m er ic al si m u la ti o n fo r th e co n ce n tr at io n of ox y ge n B (χ ) w it h an al y ti ca l so lu ti on s ob ta in ed u si n g T S M an d A D M fo r fi x ed p ar am et er va lu es α = 0. 5, β = 0. 6, γ = 0. 7, σ = 0. 5, ϕ 2 1 = 0. 5, ϕ 2 3 = 0. 7 an d va ry in g th e ϕ 2 2 ϕ 2 2 = 0 .1 ϕ 2 2 = 0 .5 ϕ 2 2 = 1 χ N u m er ic a l T S M er ro r % A D M er ro r% N u m er ic a l T S M er ro r % A D M er ro r% N u m er ic a l T S M er ro r% A D M er ro r% 0 .0 1 .0 0 1 0 0 1 .0 0 1 0 0 0 .0 0 0 0 0 1 .0 0 1 0 0 0 .0 0 0 0 0 1 .0 2 4 0 0 1 .0 2 4 0 0 0 .0 0 0 0 0 1 .0 2 4 0 0 0 .0 0 0 0 0 0 .9 7 5 9 0 0 .9 7 5 7 0 0 .0 2 0 4 9 0 .9 7 5 6 0 0 .0 3 0 7 4 0 .2 1 .0 0 0 9 6 1 .0 0 0 9 6 0 .0 0 0 0 1 1 .0 0 0 9 6 0 .0 0 0 0 1 1 .0 2 3 0 2 1 .0 2 3 0 3 0 .0 0 0 6 6 1 .0 2 3 0 3 0 .0 0 0 0 9 0 .9 7 5 9 0 0 .9 7 5 7 0 0 .0 2 0 4 9 0 .9 7 5 6 0 0 .0 3 0 7 4 0 .4 1 .0 0 0 8 4 1 .0 0 0 8 4 0 .0 0 0 0 5 1 .0 0 0 8 4 0 .0 0 0 0 3 1 .0 2 0 1 0 1 .0 2 0 1 2 0 .0 0 2 5 6 1 .0 2 0 1 0 0 .0 0 0 5 2 0 .9 7 5 9 0 0 .9 7 5 7 0 0 .0 2 0 4 9 0 .9 7 5 6 0 0 .0 3 0 7 4 0 .6 1 .0 0 0 6 5 1 .0 0 0 6 5 0 .0 0 0 1 1 1 .0 0 0 6 5 0 .0 0 0 0 3 1 .0 1 5 2 2 1 .0 1 5 2 7 0 .0 0 5 4 4 1 .0 1 5 2 4 0 .0 0 1 7 8 0 .9 7 5 9 0 0 .9 7 5 7 0 0 .0 2 0 4 9 0 .9 7 5 6 0 0 .0 3 0 7 4 0 .8 1 .0 0 0 3 8 1 .0 0 0 3 8 0 .0 0 0 1 8 1 .0 0 0 3 8 0 .0 0 0 0 2 1 .0 0 8 3 8 1 .0 0 8 4 7 0 .0 0 8 8 5 1 .0 0 8 4 3 0 .0 0 4 7 0 0 .9 7 5 9 0 0 .9 7 5 7 0 0 .0 2 0 4 9 0 .9 7 5 6 0 0 .0 3 0 7 4 1 .0 1 .0 0 0 0 3 1 .0 0 0 0 3 0 .0 0 0 2 5 1 .0 0 0 0 3 0 .0 0 0 1 9 0 .9 9 9 5 8 0 .9 9 9 7 0 0 .0 1 2 1 4 0 .9 9 9 6 8 0 .0 1 0 4 4 0 .9 7 5 9 0 0 .9 7 5 7 0 0 .0 2 0 4 9 0 .9 7 5 6 0 0 .0 3 0 7 4 0 .0 0 0 1 0 0 .0 0 0 0 3 0 .0 0 4 9 4 0 .0 0 2 9 2 0 .0 2 0 4 9 0 .0 3 0 7 4 S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 17 of 26 T a b le 3 : C om p ar is o n o f n u m er ic a l si m u la ti on fo r th e co n ce n tr at io n of p ro d u ct C (χ ) w it h an al y ti ca l so lu ti on s ob ta in ed u si n g T S M an d A D M fo r fi x ed p ar am et er va lu es α = 0. 5, β = 0. 6, γ = 0. 7, σ = 0. 5, ϕ 2 1 = 0. 5, ϕ 2 2 = 0. 6 an d va ry in g th e ϕ 2 3 ϕ 2 3 = 0 .3 ϕ 2 3 = 0 .5 ϕ 2 3 = 0 .7 χ N u m er ic a l T S M er ro r % A D M er ro r% N u m er ic a l T S M er ro r % A D M er ro r% N u m er ic a l T S M er ro r% A D M er ro r% 0 .0 0 .9 9 5 6 0 0 .9 9 5 6 0 0 .0 0 0 0 0 1 .0 0 0 0 0 0 .0 0 0 0 0 0 .9 8 7 7 0 0 .9 8 7 8 0 0 .0 1 0 1 2 0 .9 8 7 8 0 0 .0 1 0 1 2 1 .0 2 4 0 0 1 .0 2 4 0 0 0 .0 0 0 0 0 1 .0 2 4 0 0 0 .0 0 0 0 0 0 .2 0 .9 9 5 7 8 0 .9 9 5 7 8 0 .0 0 0 2 2 1 .0 0 0 0 0 0 .0 0 0 0 0 0 .9 8 8 1 9 0 .9 8 8 2 9 0 .0 0 9 4 6 0 .9 8 8 2 9 -0 .0 1 0 1 7 1 .0 2 3 0 2 1 .0 2 3 0 3 0 .0 0 0 6 6 1 .0 2 3 0 3 0 .0 0 0 0 9 0 .4 0 .9 9 6 3 1 0 .9 9 6 3 0 0 .0 0 0 8 9 1 .0 0 0 0 1 0 .0 0 0 0 0 0 .9 8 9 6 7 0 .9 8 9 7 4 0 .0 0 7 4 0 0 .9 8 9 7 7 0 .0 1 0 3 6 1 .0 2 0 1 0 1 .0 2 0 1 2 0 .0 0 2 5 6 1 .0 2 0 1 0 0 .0 0 0 5 2 0 .6 0 .9 9 7 1 9 0 .9 9 7 1 8 0 .0 0 1 9 3 1 .0 0 0 0 1 0 .0 0 0 0 0 0 .9 9 2 1 3 0 .9 9 2 1 7 0 .0 0 4 2 7 0 .9 9 2 2 4 0 .0 1 0 8 9 1 .0 1 5 2 2 1 .0 1 5 2 7 -0 .0 0 5 4 4 1 .0 1 5 2 4 0 .0 0 1 7 8 0 .8 0 .9 9 8 4 3 0 .9 9 8 4 0 0 .0 0 3 1 4 1 .0 0 0 0 2 0 .0 0 0 0 0 0 .9 9 5 5 7 0 .9 9 5 5 8 0 .0 0 0 6 1 0 .9 9 5 6 9 0 .0 1 2 0 6 1 .0 0 8 3 8 1 .0 0 8 4 7 0 .0 0 8 8 5 1 .0 0 8 4 3 0 .0 0 4 7 0 1 .0 1 .0 0 0 0 3 0 .9 9 9 9 8 0 .0 0 4 2 9 1 .0 0 0 0 4 0 .0 0 0 0 0 1 .0 0 0 0 0 0 .9 9 9 9 7 0 .0 0 2 8 2 1 .0 0 0 1 4 0 .0 1 4 3 5 0 .9 9 9 5 8 0 .9 9 9 7 0 0 .0 1 2 1 4 0 .9 9 9 6 8 -0 .0 1 0 4 4 0 .0 0 1 7 4 0 .0 0 0 0 0 0 .0 0 4 8 4 0 .0 1 1 3 3 0 .0 0 4 9 4 0 .0 0 2 9 2 S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 18 of 26 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless Distance 0.975 0.98 0.985 0.99 0.995 1 D im en si o n le ss C o n ce n tr a ti o n o f O x y g en =0.5 = 10 = 100 (a) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Dimensionless Distance 0.98 0.985 0.99 0.995 1 1.005 D im en si o n le ss C o n ce n tr a a a ti o n o f O x y g en = 10 = 100 = 0.5 (b) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless Distance 0.975 0.98 0.985 0.99 0.995 1 1.005 D im en si o n le ss c o n ce n tr a ti o n o f O x y g en = 10 = 0.7 = 100 (c) 0 0.2 0.4 0.6 0.8 1 Dimensionless Distance 0.975 0.98 0.985 0.99 0.995 1 1.005 D im en si o n le ss c o n ce n tr a ti o n o f O x y g en = 100 = 10 = 0.5 (d) Figure 4: Dimensionless concentration for Oxygen B(χ) versus the dimensionless distance χ for various values of parameters (4a) α, (4b) β, (4c) γ, (4d) σ where ‘· · · ’ denotes the numerical solution, ‘***’ represents the ADM solution, and ‘⋄⋄⋄’ corresponds to the TSM solution. S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 19 of 26 0 10 20 30 40 50 1.82 1.84 1.86 1.88 1.9 1.92 1.94 1.96 1.98 2 E ff e c ti v e n e s s F a c to r 2 2 = 1 2 2 = 0.8 2 2 = 0.6 2 2 = 0.4 2 2 = 0.2 (a) 0 10 20 30 40 50 1.85 1.9 1.95 2 E ff e c ti v e n e s s F a c to r 2 2 = 0.8 2 2 = 0.6 2 2 = 0.4 2 2 = 0.2 2 2 = 1 (b) 0 10 20 30 40 50 1.8 1.82 1.84 1.86 1.88 1.9 1.92 1.94 1.96 1.98 2 E ff e c ti v e n e s s F a c to r 2 2 = 1 2 2 = 0.8 2 2 = 0.6 2 2 = 0.4 2 2 = 0.2 (c) 0 5 10 15 20 25 30 35 40 45 50 1.85 1.9 1.95 2 E ff e c ti v e n e s s F a c to r 2 2 = 0.2 2 2 = 0.4 2 2 = 0.6 2 2 = 0.8 2 2 = 1 (d) Figure 5: Effectiveness factor versus saturation parameters (5a) σ, (5b) α, (5c) β, (5d) γ for various values of diffusion parameter ϕ2 2 6. Conclusion The mathematical model of the carbohydrate oxidation by oxygen catalyzed by the bienzyme glucose dehydrogenase/laccase system immobilized into MR with carbon nan- otubes is analyzed by utilizing the semi-analytical methods, namely TSM and ADM. When compared with the numerical simulation, satisfactory results are obtained. The obtained semi-analytical solutions are used to analyze the effects of diffusion and saturation pa- rameters on the concentrations and effective working of the batch reactor. It is observed that the diffusion parameter ϕ2 2 has the major effect on the working of the batch reactor. As the increase in the diffusion parameter increases, the effectiveness factor of the batch reactor increases. S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 20 of 26 The proposed theoretical model for the oxidation of carbohydrates by oxygen, catalyzed by the bienzyme glucose dehydrogenase/lactase system immobilized within an MR incor- porating carbon nanotubes, offers significant utility for experimental researchers. This model aids in optimizing batch reactor design by evaluating enzyme catalytic loading and particle size, providing a valuable tool for analyzing reaction conditions. To extend the present mathematical model of carbohydrate oxidation process catalyzed by the glucose dehydrogenase/laccase bienzyme system within a carbon nanotube-based micro-reactor, future work can focus on the development of time-dependent non-steady- state model of the micro-reactor system. This model would enable a deeper understanding of the transient dynamics associated with substrate consumption and product formation, especially during initial stages or under varying input conditions. Further more this model can be extended to the multi-dimensional spatial models for more accurate representation of the porous micro-environment with spherical micro-bioreactors. A. Semi-analytical expressions utilizing Taylor Series Method The solution for the governing, Eqs. (12) - (14), for concentration using TSM is ex- pressed as follows: Considering the series solution for lactose around χ = 0, we have: A(χ) = ∞∑ i=0 diA dχi ∣∣∣∣∣ χ=1 (26) = A(1) + χ− 1 1! dA dχ ∣∣∣∣ χ=1 + (χ− 1)2 2! d2A dχ2 ∣∣∣∣ χ=1 + (χ− 1)3 3! d3A dχ3 ∣∣∣∣ χ=1 + · · · (27) considering the substrate at χ = 1 A(χ)|χ=1 = 1 (28) dA dχ ∣∣∣∣ χ=1 = l (29) d2A dχ2 ∣∣∣∣ χ=1 = −2l + 2ϕ2 1 α+ 2β + γ + 2σ (30) d3A dχ3 ∣∣∣∣ χ=1 = (( 6α2 + (24β + 12γ + 24σ)α+ 24β2 + (4ϕ2 1 + 24γ + 48σ)β+ 6(γ + 2σ)2 ) l + 2ϕ2 1 ((n− 2)α− 4β − 2γ − 4σ) ) (α+ 2β + γ + 2σ)2 (31) S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 21 of 26 d4A dχ4 ∣∣∣∣ χ=1 =  − 8(γ + 2σ + α)βϕ2 1l 2 + ( 16β ((n− 1)α− 2β − γ − 2σ)ϕ2 1 −192 (γ 2 + σ + α 2 + β )3) l − 8ϕ2 1 ( −βϕ2 1 + (n− 2)α2+( (n2 + 2n− 8)β + ( n2 2 + n− 4 ) γ + (n2 + 2n− 8)σ − ϕ2 3 4 ) α −8 (γ 2 + σ + β )2)  (α+ 2β + γ + 2σ)3 (32) d5A dχ5 ∣∣∣∣ χ=1 =  16β ( (−3l + 7n 2 − 2)α+ (l − 4)β − 3 ( l + 2 3 ) (γ + 2σ) ) ϕ4 1 + ( (40n− 80)α3 + ( (24l3 + (−48n+ 64)l2 + (24n2 − 128n+ 80)l +64n2 + 160n− 480)β + (32n2 + 80n− 240)γ + (64n2 + 160n− 480)σ +2ϕ2 3(n− 4) ) α2 + ( ((48n+ 128)l2 + (−96n2 − 256n+ 320)l + 48n3 +128n2 + 160n− 960)β2 + ( (48γ + 96σ)l3 − 48(n− 8 3 )(γ + 2σ)l2 + ( (−48n2 − 128n+ 160)γ + (−96n2 − 256n+ 320)σ + 28ϕ2 3 ) l +(48n3 + 128n2 + 160n− 960)γ + (96n3 + 256n2 + 320n− 1920)σ −24ϕ2 3(n+ 2 3 ) ) β + 12 ( (n3 + 8 3 n2 + 10 3 n− 20)γ + ( 2n3 + 16 3 n2 + 20 3 n −40)σ − ϕ2 3(n+ 2 3 ) ) (γ + 2σ) ) α+ (320l − 640)β3 + 128 ( l2 + 5 2 l − 15 2 ) (γ + 2σ)β2 + 24(l3 + 8 3 l2 + 10 3 l − 20)(γ + 2σ)2β − 80(γ + 2σ)3 ) ϕ2 1 + 1920 (γ 2 + σ + α 2 + β )4 l  (α+ 2β + γ + 2σ)4 (33) using the Eqs. (28) - (33) we get the solution of the lactose concentration as given in (18). By following the same methodology we can find the semi-analytical solutions of product B(χ) and oxygen C(χ) concentrations. B. Semi-analytical expressions utilizing Adomian Decomposition Method The non linear governing Eqs. (12) - (14) are solved using the ADM, by applying L[A(χ)] = ϕ2 1N [A(χ), C(χ), B(χ)], (34) L[C(χ)] = −ϕ2 2N [A(χ), C(χ), B(χ)], (35) L[B(χ)] = ϕ2 3 2 N [A(χ), C(χ), B(χ)], (36) S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 22 of 26 where L = d2 dχ2 , L−1(.) = ∫ ∫ (.) dχdχ, N [A(χ)] = − 2 χ dA(χ) dχ = ϕ2 1 2A(χ)B(χ) αA(χ) + 2βB(χ) + γA(χ)B(χ) + 2σA(χ)B(χ) , (37) utilizing the inverse operator in Eqs. (34) - (36) A(χ) = D1χ+D2 + ϕ2 1L −1N [A(χ)], (38) C(χ) = D3χ+D4 + ϕ2 2L −1N [C(χ)], (39) B(χ) = D5χ+D6 + ϕ2 3L −1N [B(χ)], (40) where D1, D2, D3, D4, D5, and D6 are integration constants. The ADM assumes an infinite series solution for the unknown functions A(x), B(x), C(x) of Eqs. (34)–(36), given by [32–34] as follows: A(x) = ∞∑ n=0 An(x), (41) C(x) = ∞∑ n=0 Cn(x), (42) B(x) = ∞∑ n=0 Bn(x), (43) and decomposing the non-linear operator N as N [A(x)] = ∞∑ n=0 Fn(x), (44) Fn are called Adomian polynomials of A0, A1, . . . , An [33, 34] and are formed by Fn(x) = 1 n! dn dλn [ N ( ∞∑ i=0 λiAi(x) )]∣∣∣∣∣ λ=0 , n = 0, 1, 2, . . . (45) Substituting Eqs. (41)–(45) into Eqs. (38)–(40), we get the following iterations: A1 = 1 A′ 0 = 0 (46) An+1 = D2 + ϕ2 1L −1N [A(χ)] = D2 + ϕ2 1L −1Fn(χ), n ≥ 0 (47) Similarly for, Cn+1 = D4 + ϕ2 2L −1N [B(χ)] = D4 + ϕ2 2L −1Fn(χ), n ≥ 0 (48) S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 23 of 26 Bn+1 = D6 + ϕ2 3L −1N [C(χ)] = D6 + ϕ2 3L −1Fn(χ), n ≥ 0 (49) F0(χ) = N [A(χ)] = A0B0 αA0 + 2βB0 + γA0B0 + σA0B0 = 1 3(α+ 2β + γ + 2σ) , (50) using Eqs. (47) - (49) we get A1(χ) = D2 + ϕ2 1L −1F0(χ) = D2 + ϕ2 1L −1 ( 1 3(α+ 2β + γ + 2σ) ) , = ϕ2 1(χ 2 − 1) 3(α+ 2β + γ + 2σ) (51) Similarly for C1(χ) = D4 + ϕ2 2L −1F0(χ) = D4 + ϕ2 2L −1 ( 1 3(α+ 2β + γ + 2σ) ) , = ϕ2 2(1− χ2) 3(α+ 2β + γ + 2σ) , (52) B1(χ) = D6 + ϕ2 3L −1F0(χ) = D6 + ϕ2 3L −1 ( 1 3(α+ 2β + γ + 2σ) ) , = ϕ2 3(χ 2 − 1) 6(α+ 2β + γ + 2σ) , (53) Then F1(χ) can be obtained by following boundary conditions A1(0) = 0, A1(1) = 1 (54) F1(χ) = d dλ (N [A0 + λA1])λ=0 and A′ i(0) = 0, Ai(1) = 0, i ≥ 0 (55) A2(χ) = D2 + ϕ2 1L −1F0(χ), = D2 + ϕ2 1L −1 ( (χ2 − 1)(αϕ2 3 + 4βϕ2 1) , 6(α+ γ + 2β + 2σ)3 ) A2(χ) = ϕ2 1 (χ2 − 1)(αϕ2 3 + 4βϕ2 1) 6(α+ γ + 2β + 2σ)3 ( χ4 20 − χ2 6 + 7 60 ) , (56) the concentration of lactose A(χ) is determined by summing Eqs. (46), (53), and (56). Using a similar approach, the solutions for the concentration of the product C(χ) and the concentration of oxygen B(χ) can also be obtained. Acknowledgements The authors are very much thankful to the management, SRM Institute of Science and Technology for their continuous support and encouragement. S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 24 of 26 Author Contributions Daniel S; Data curation, Writing-original draft, Formal analysis, Methodology, Soft- ware, Visualization. Mallikarjuna M; Data curation, Writing-original draft, Formal anal- ysis, Methodology, Software, Visualization. Senthamarai R; Formal analysis, Method- ology, Validation, Resources, Writing review and editing and Supervision. Declarations Funding This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. Conflict of interest The authors declare that they have no conflict of interest. References [1] J. Iqbal, S. Iqbal, and C. E. Müller. Advances in immobilized enzyme microbioreactors in capillary electrophoresis. Analyst, 138(11):3104–3116, 2013. [2] D. Schäpper, M. N. H. Z Alam, N. Szita, A. Eliasson Lantz, and K. V Gernaey. Appli- cation of microbioreactors in fermentation process development: a review. Analytical and bioanalytical chemistry, 395:679–695, 2009. [3] G. Maria. Enzymatic reactor selection and derivation of the optimal operation pol- icy, by using a model-based modular simulation platform. Computers & chemical engineering, 36:325–341, 2012. [4] E. Papadakis, S. Pedersen, A. K. Tula, M. Fedorova, M. J. Woodley, and R. Gani. Model-based design and analysis of glucose isomerization process operation. Com- puters & Chemical Engineering, 98:128–142, 2017. [5] W. Tischer and F. Wedekind. Immobilized enzymes: methods and applications. Biocatalysis-from discovery to application, pages 95–126, 1999. [6] A. Y. Khan, S. B. Noronha, and R. Bandyopadhyaya. Glucose oxidase enzyme im- mobilized porous silica for improved performance of a glucose biosensor. Biochemical engineering journal, 91:78–85, 2014. [7] L. T. Zhu, W. Y. Ma, and Z. H. Luo. Influence of distributed pore size and porosity on mto catalyst particle performance: Modeling and simulation. Chemical Engineering Research and Design, 137:141–153, 2018. [8] P. M. Doran. Bioprocess engineering principles. Elsevier, 1995. [9] J. Villadsen. Principles of metabolic engineering. Biotechnology, 3:226–256, 2009. [10] M. Al-Shannag, Z. Al-Qodah, J. Herrero, A. J. Humphrey, and F. Giralt. Using a wall-driven flow to reduce the external mass-transfer resistance of a bio-reaction system. Biochemical engineering journal, 39(3):554–565, 2008. [11] R. Baronas, J. Kulys, and L. Petkevičius. Modeling carbohydrates oxidation by oxy- gen catalyzed by bienzyme glucose dehydrogenase/laccase system immobilized into S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 25 of 26 microreactor with carbon nanotubes. Journal of Mathematical Chemistry, 59:168– 185, 2021. [12] J. B. Sánchez, M. Vuono, and D. Dionisi. Model-based comparison of sequencing batch reactors and continuous-flow activated sludge processes for biological wastew- ater treatment. Computers & Chemical Engineering, 144:107127, 2021. [13] F. F. Rivera, T. Pérez, L. F. Castañeda, and J. L. Nava. Mathematical modeling and simulation of electrochemical reactors: A critical review. Chemical engineering science, 239:116622, 2021. [14] O. Harker-Sanchez, A. A. Jaramillo, and D. M. Arias. Method to obtain parameters k2, k3 for dilution rate observer in am2 model of the anaerobic digestion process in a batch reactor. Energy Sources, Part A: Recovery, Utilization, and Environmental Effects, 46(1):3110–3123, 2024. [15] M. Mallikarjuna and R. Senthamarai. An amperometric biosensor and its steady state current in the case of substrate and product inhibition: Taylors series method and adomian decomposition method. Journal of Electroanalytical Chemistry, 946:117699, 2023. [16] M. Mallikarjuna, J. L. G. Guiraob, and R. Senthamarai. Theoretical analysis of am- perometric biosensor with substrate and product inhibition involving non-michaelis- menten kinetics. MATCH Communications in Mathematical and Computer Chem- istry, 93:319–347, 2025. [17] M. Mallikarjuna and R. Senthamarai. Mathematical analysis of batch reactor per- formance for the enzymatic synthesis of cephalexin: Laplace homotopy perturbation method and adomian decomposition method. Partial Differential Equations in Ap- plied Mathematics, 11:100806, 2024. [18] M. Sivakumar, M. Mallikarjuna, and R. Senthamarai. A kinetic non-steady-state analysis of immobilized enzyme systems without external mass transfer resistance. International Journal of Analysis and Applications, 22:31–31, 2024. [19] M. Sivakumar, M. Mallikarjuna, and R. Senthamarai. A kinetic non-steady state analysis of immobilized enzyme systems with external mass transfer resistance. AIMS Mathematics, 9(7):18083–18102, 2024. [20] W. Al-Hayani. Combined laplace transform-homotopy perturbation method for sine- gordon equation. Applied Mathematics & Information Sciences, 10:1781–1786, 2016. [21] R. Senthamarai and T. N. Saibavani. Substrate mass transfer: analytical approach for immobilized enzyme reactions. In Journal of Physics: Conference Series, volume 1000, page 012146. IOP Publishing, 2018. [22] W. M. Al-Hayani and M. T. Younis. Solving fuzzy system of boundary value problems by homotopy perturbation method with green’s function. European Journal of Pure and Applied Mathematics, 16(2):1236–1259, 2023. [23] W. M. Al-Hayani and M. T. Younis. The homotopy perturbation method for solving nonlocal initial-boundary value problems for parabolic and hyperbolic partial differen- tial equations. European Journal of Pure and Applied Mathematics, 16(3):1552–1567, 2023. [24] A. El-Ajou, H. Al-ghananeem, R. Saadeh, A. Qazza, and M. A. N. Oqielat. A modern S. Daniel, M. Mallikarjuna, R. Senthamarai / Eur. J. Pure Appl. Math, 18 (3) (2025), 6203 26 of 26 analytic method to solve singular and non-singular linear and non-linear differential equations. Frontiers in Physics, 11:1167797, 2023. [25] M. O. Olayiwola, A. I. Alaje, A. O. Yunus, K. A. Adedokun, and K. A. Bashiru. A mathematical modeling of covid-19 treatment strategies utilizing the laplace adomian decomposition method. Results in Control and Optimization, 14:100384, 2024. [26] R. O. Awonusika. Analytical solution of a class of lane–emden equations: Adomian decomposition method. The Journal of Analysis, 32(2):1009–1056, 2024. [27] A. S. Mahmood, L. Casasús, and W. Al-Hayani. Analysis of resonant oscillators with the adomian decomposition method. Physics Letters A, 357(4-5):306–313, 2006. [28] W. Al-Hayani. Adomian decomposition method with green’s function for solving twelfth-order boundary value problems. Applied Mathematical Sciences, 9:353–368, 01 2015. [29] M. Usman, H. U. Khan, M. Sarwar, N. Fatima, and K. Abodayeh. An innovative analytical result of two-dimensional heat equation using joint mechanism of natural transform and adomian decomposition method. European Journal of Pure and Applied Mathematics, 18(2):6051–6051, 2025. [30] V. Silambuselvi, P. Jeyabarathi, N. Jha, K. Angaleeswari, T. R. K. Kumar, and L. Ra- jendran. Theoretical analysis of facilitated diffusion process in a liquid membrane: Adomian decomposition method. International Journal of Electrochemical Science, 19(12):100855, 2024. [31] B. Mebarek, A. Maatoug, S. A. M. Mostefaoui, H. Benali, and Y. El Guerri. Ado- mian decomposition method for modelling the growthof feb/fe2b layer in boronizing process. Zastita Materijala, 66(1):148–157, 2025. [32] G. Adomian. Nonlinear stochastic operator equations. Academic press, 1986. [33] G. Adomian. Nonlinear stochastic systems theory and applications to physics, vol- ume 46. Springer Science & Business Media, 1988. [34] G. Adomian. Solving frontier problems of physics: the decomposition method. Springer Science & Business Media, 1994.