Numerical Modeling and Investigation of Amperometric Biosensors with Perforated Membranes

The present paper aims to investigate the influence of perforated membrane geometry on the performance of biosensors. For this purpose, a 2-D axisymmetric model of an amperometric biosensor is analyzed. The governing equations describing the reaction-diffusion equations containing a nonlinear term related to the Michaelis–Menten kinetics of the enzymatic reaction are introduced. The partial differential governing equations, along with the boundary conditions, are first non-dimensionalized by using appropriate dimensionless variables and then solved in a non-uniform unstructured grid by employing the Galerkin Finite Element Method. To examine the impact of the hole-geometry of the perforated membrane, seven different geometries—including cylindrical, upward circular cone, downward circular cone, upward paraboloid, downward paraboloid, upward concave paraboloid, and downward concave paraboloid—are studied. Moreover, the effects of the perforation level of the perforated membrane, the filling level of the enzyme on the transient and steady-state current of the biosensor, and the half-time response are presented. The results of the simulations show that the transient and steady-state current of the biosensor are affected by the geometry dramatically. Thus, the sensitivity of the biosensor can be influenced by different hole-geometries. The minimum and maximum output current can be obtained from the cylindrical and upward concave paraboloid holes. On the other hand, the least half-time response of the biosensor can be obtained in the cylindrical geometry.


Introduction
The detection of bio-particles (like nucleic acids, proteins, bacteria, and viruses) plays a crucial role in the early diagnosis and prevention of various diseases, such as breast and prostate cancers. Nowadays, the start of diseases and their progress can be rapidly assessed using biosensors [1,2]. A quick and precise detection can lead to the administration of appropriate medicine and reduces the possibility of reaching a critical stage of the disease [3]. Therefore, the study of such biosensors has simulation and revealed that the permeability to the analyte of the external membrane and its thickness are the most influential parameters for enhancing the response time of sandwich-type biosensors.
An outer porous or perforated membrane is mainly utilized in practical biosensors for increasing their stability and also prolonging their calibration curves [16][17][18]. Schulmeister and Pfeiffer [19] developed the mathematical modeling of a biosensor in the presence of a perforated membrane. However, the geometry of the membrane was not considered in their 1-D model. The two-dimensional-in-space mathematical modeling of biosensors was proposed by Baronas et al. [20,21]. In [21], Baronas et al. investigated the influence of the presence of a cylindrical perforated membrane on the performance of the biosensor. They assumed that the perforated membrane was completely filled with enzymes and showed that the perforation level of the membrane could strongly affect the steady-state and also half-time response of the biosensor. In the other study, Baronas et al. [20] surveyed the influence of the filling level of enzymes in the perforated membrane and revealed that the impact of the filling level on the biosensor response diminished as the radius of the hole decreased.
It is worth noting that in comparison to 1-D cases, the 2-D modeling of biosensors is highly computationally demanding and may not be efficient in some situations. As such, in another study, Baronas et al. [22] evaluated the conditions in which 1-D models are accurate enough to be employed for the mathematical modeling of biosensors. They found that 1-D models cannot be used when the hole of the perforated membrane is very small; thus, in such cases, 2-D models should be utilized to obtain accurate outcomes. In addition to this, the geometry of the perforated membrane governs the mass transport by the diffusion process and thus can affect the overall performance of amperometric biosensors. Previous surveys have significantly increased knowledge about the impact of influential parameters on the performance of biosensors. However, in all of the studies to date, only cylindrical perforated membranes have been analyzed and, to the best of the authors' knowledge, the impact of the shape of the mentioned membrane on the steady-state current and time response of biosensors has not been addressed yet.
With the above literature review, it is now clear that this study aims to numerically survey the influence of the geometry of the perforated membrane on the performance of biosensors. For this purpose, a 2-D axisymmetric amperometric biosensor with selective and perforated membranes was analyzed. Seven different geometries were designed and modeled. The governing equations were solved using the Galerkin Finite Element Method. The results were first compared and validated with four different works. Then, the influence of the geometry of the perforated membranes on the steady-state current and also the time response of the biosensor was discussed.

Structure of the Biosensor
The schematic view of the biosensor is provided in Figure 1. It is assumed that the thickness of the perforated and selective membranes is much lower than the length and width of the biosensor. The holes in the perforated membrane can be modeled to form the same unit cells with hexagonal patterns. Only half of the cross-section of the mentioned unit cell is studied because of its symmetry [20,21]. Scheme 1 illustrates various regions of the biosensor. While in the bulk solution region both convection and diffusion mechanisms (without reaction) can be found, the diffusion mode governs the other layers [14,23]. The enzyme-catalyzed reaction takes place in the enzyme region. In the selective membrane, the production detection process is performed. Finally, in the electrode, the biochemical energy from the enzyme reaction is converted by the transducer to the output current. Scheme 1 illustrates various regions of the biosensor. While in the bulk solution region both convection and diffusion mechanisms (without reaction) can be found, the diffusion mode governs the other layers [14,23]. The enzyme-catalyzed reaction takes place in the enzyme region. In the selective membrane, the production detection process is performed. Finally, in the electrode, the biochemical energy from the enzyme reaction is converted by the transducer to the output current.

Formulation of the Problem
The biosensor can be mathematically modeled in a 2-D domain consisting of the discussed regions. The Michaelis-Menten model is utilized for the study of the kinetics of enzymes [24]. According to the Michaelis-Menten model, the substrate (S) binds to the enzyme (E) with the reaction rate of k1 and converts to the ES. The formed complex (ES) then dissociates during the second reaction with the reaction rate of k2 and produces the product (P), regenerating the enzyme [20]. It is worth noting that the rate of forward reaction rate (k2) for 2 k ES E P    is much higher than the reverse one, and thus the reverse reaction can be neglected. Following Baronas et al. [25], the reaction model can be shown as follows: Since the concentration of the substrate is much more than the ES, it is reasonable to consider a constant concentration for the ES. Thus, the concentration of the enzyme production can be calculated using the following equation [26]:  Scheme 1 illustrates various regions of the biosensor. While in the bulk solution region both convection and diffusion mechanisms (without reaction) can be found, the diffusion mode governs the other layers [14,23]. The enzyme-catalyzed reaction takes place in the enzyme region. In the selective membrane, the production detection process is performed. Finally, in the electrode, the biochemical energy from the enzyme reaction is converted by the transducer to the output current.

Formulation of the Problem
The biosensor can be mathematically modeled in a 2-D domain consisting of the discussed regions. The Michaelis-Menten model is utilized for the study of the kinetics of enzymes [24]. According to the Michaelis-Menten model, the substrate (S) binds to the enzyme (E) with the reaction rate of k1 and converts to the ES. The formed complex (ES) then dissociates during the second reaction with the reaction rate of k2 and produces the product (P), regenerating the enzyme [20]. It is worth noting that the rate of forward reaction rate (k2) for is much higher than the reverse one, and thus the reverse reaction can be neglected. Following Baronas et al. [25], the reaction model can be shown as follows: Since the concentration of the substrate is much more than the ES, it is reasonable to consider a constant concentration for the ES. Thus, the concentration of the enzyme production can be calculated using the following equation [26]:

Formulation of the Problem
The biosensor can be mathematically modeled in a 2-D domain consisting of the discussed regions. The Michaelis-Menten model is utilized for the study of the kinetics of enzymes [24]. According to the Michaelis-Menten model, the substrate (S) binds to the enzyme (E) with the reaction rate of k 1 and converts to the ES. The formed complex (ES) then dissociates during the second reaction with the reaction rate of k 2 and produces the product (P), regenerating the enzyme [20]. It is worth noting that the rate of forward reaction rate (k 2 ) for ES k 2 → E + P is much higher than the reverse one, and thus the reverse reaction can be neglected. Following Baronas et al. [25], the reaction model can be shown as follows: Since the concentration of the substrate is much more than the ES, it is reasonable to consider a constant concentration for the ES. Thus, the concentration of the enzyme production can be calculated using the following equation [26]: where S, E, and E S represent the concentrations of the substrate, enzyme, and the enzyme-substrate complex, respectively. E 0 is the initial concentration of the enzyme. By replacing the expression E = E 0 − E S in Equation (2), the production rate (v) can be written as follows [25,27]: In Equation (3), P represents the concentration of the product. Moreover, the reaction rate approaches the saturated reaction rate when the E 0 = E S . For this case, V m = k 2 E 0 represents the maximum reaction rate and K m = (k 2 + k −1 )/k 1 is the Michaelis-Menten constant, showing the concentration of the substrate when the reaction velocity is equal to one half of the maximal velocity of the reaction. The mass diffusion in the layers can be expressed over time with Fick's law of diffusion in two or three dimensions, as follows: where t is the time, ∆ is the Laplace operator, c is the concentration of a component, and D is the corresponding diffusion coefficient. The value of the diffusion coefficient indicates the quality of the transfer, which means that the larger the diffusion coefficient, the better the diffusion mass transfer. Shown in Figure 2 are the seven studied unit cells with different geometrical designs. Except for the cylindrical one (Figure 2a), the perforated membranes were designed to be in upward or downward pairs of a circular cone, paraboloid, and concave paraboloid. As such, the minimum and maximum radii of the perforated membranes are set to be a * 2 and a * 3 , respectively. Moreover, a * 1 corresponds to the hole radius. The thickness of the selective membrane is b * are the perforated membrane thickness and the diffusion layer thickness, respectively. The filling level of enzymes is shown with b 3 . Region Ω 1 indicates the selective membrane, Ω 2 specifies the layer occupied with the enzymes, Ω 3 is the external diffusion layer, and Ω 4 is the impermeable carrier of the perforated membrane.
Sensors 2020, 20, x 6 of 17 It is assumed that the selected layer has a uniform thickness, and the thicknesses of both the perforated membrane and the selective membrane in the biosensor are much less than its length and width. The volume of the bulk solution is considered to be extremely high, while the substrate concentration can be assumed to be constant [28,29]. The concentration of the ES form can be assumed to be constant because the biosensor reached the steady-state condition in a short time. By using Equation (2), the maximum rate of enzymatic reaction can be obtained. Therefore, according to Figure   It is assumed that the selected layer has a uniform thickness, and the thicknesses of both the perforated membrane and the selective membrane in the biosensor are much less than its length and width. The volume of the bulk solution is considered to be extremely high, while the substrate concentration can be assumed to be constant [28,29]. The concentration of the ES form can be assumed to be constant because the biosensor reached the steady-state condition in a short time. By using Equation (2), the maximum rate of enzymatic reaction can be obtained. Therefore, according to Figure 2, the biosensor described from the process starts (t* > 0) with the system of reaction-diffusion equations as follows: are the product and substrate concentrations, respectively. The symbols r* and z* indicate the radial and axial distances, respectively, and t* denotes elapsed time. P * 1 is the product concentration on the surface of the electrode; this concentration and its variation are important for measuring the change in the current of the biosensor (see Equation (26)). Here, P * 2 and P * 3 represent the product concentrations of the enzyme and external diffusion layers, respectively. Since the reaction product oxidizes at the electrode surface rapidly, it can be assumed that its concentration at the electrode surface (z * = 0) is zero [21].
At the interface between the bulk solution and the biosensor (z * = b * 5 ), the concentration of the enzyme reaction is assumed to be zero because the reaction has not yet occurred. S * 2 and S * 3 are the substrate concentration on the selective membrane and the perforated membrane. D * 1 , D * 2 , and D * 3 represent the diffusion coefficient in these regions. The initial conditions are described as (t * > 0): where S * 0 is the substrate concentration in the solution, Ω i is the area of each layer, and Γ i is the upper surface of Ω i . Due to the oxidation of the reaction product at the electrode surface, its concentration at the electrode surface is considered to be zero (t * > 0): Since no reaction has occurred on the boundary, (z * = b * 5 )Γ 3 , the concentration of the enzymatic product, has to be constant and zero. During the biosensor operation, the substrate concentration, as well as the product over the enzyme surface (bulk solution and membrane interface), remains constant; therefore (t * > 0): The non-leakage condition for the symmetric boundaries of the unit cell is considered as follows: where n denotes the normal vector of the hole surface. At the boundary between the two regions, it is necessary to apply matching conditions. According to the conservation of mass, the mass flux of particles at the interfaces is equal on both sides. In addition to this, it can be assumed that the substrate concentration and the enzyme product on the boundaries are equal for both sides: In physical experiments, the anodic current generated at the electrode surface is considered as the amperometric biosensor response, which is a linear function of the reaction product gradient at the surface of the electrode. According to the Faraday and Fick laws, the current density on the electrode surface can be obtained as: where n e is the number of electrons involved in a charge transfer at the electrode surface and F is the Faraday constant [4]. The steady-state current I* can be obtained when the concentration of the production remains constant at the sensor surface: Finally, in the present study, T 0.5 , half of the required time to reach the steady-state current, is used to evaluate the dynamics of the biosensor operation [25]: Sensors 2020, 20, 2910 8 of 16

Dimensionless Mathematical Model
The governing Equations (6) and the corresponding boundary and initial conditions (7)-(11) are non-dimensionalized by introducing the following set of dimensionless variables: Substituting Equation (16) into the governing Equation (5), the dimensionless nonlinear reaction-diffusion equations can be summarized as follows: where D i = D * i /D 1 is the dimensionless diffusion coefficient and σ 2 = V max a * 1 2 /K m D 1 is the diffusion modulus (Damköhler number), which compares the rate of enzyme reaction (V max /K m ) with the mass transport through the enzyme layer (D 1 /a * 1 2 ) [25]. In fact, for high values of σ 2 , the diffusion mechanism is considered to control the biosensor response, while, when it approaches zero, the enzyme kinetics prevail in the response [6,21]. The non-dimensional initial conditions are as follows: and the dimensionless boundary conditions read: The non-dimensional matching conditions are: and the non-dimensional density of the current is given by: and the dimensionless output current reads as: To compare the effects of hole-geometry and filling level on the i(t) and I in different geometries, the following parameters are used [20][21][22]: where γ is the level of the enzyme filling within the holes. For instance, γ = 0.5 represents a case where half of the perforated membrane is filled with the enzyme. α represents the perforation level of the perforated membrane. For the set of studied geometrical designs, the maximum radius of the holes (a 3 ) is considered to be constant, and thus α is a measure of the volume fraction of the holes in the perforated membrane. In other words, increasing the perforation level (α) augments the surface area of the impermeable part of the perforated membrane (Ω 4 ).

Numerical Approach
The partial differential governing Equation (17) and the corresponding initial and boundary conditions (18)- (20) are solved numerically utilizing the Galerkin method. Using the basis set {ξ k } N k=1 , the product and the enzyme concentrations can be expanded as follows: It should be noted that the function is the same for both variables. By substituting in Equation (17), the following residual equations can be obtained at the internal domain nodes: Sensors 2020, 20, 2910 10 of 16 In order to numerically evaluate the integral terms, the bi-quadratic function with three-point Gaussian quadrature is employed. Then, the Newton-Raphson method is used to evaluate the coefficients P jk ( j = 1, 2, 3) and S jk ( j = 1, 2) in the nonlinear residual Equation (26). More details about the solution procedure can be found in the Reddy studies [30]. The iterative process is stopped when the convergence criterion (R j i ) 2 ≤ 10 −7 1 ≤ j ≤ 5 is reached. Moreover, a time step of 0.01 s is used to evaluate the half-time response of the biosensor.

Validation of Computations
To ensure about the validity of the solution, the results of the present study have been compared against the different numerical and experimental results of previously published surveys. Shown in Figure 3 are the comparisons of the results of the present study with those provided by Baronas et al. [10] for a 1-D case without perforated and selective membranes (D s = D p = 3.010 −6 cm 3 /s, K m = 10 −7 cm 3 /s). Baronas et al. [10] employed FDM for solving the set of non-linear equations. Further, to compare the accuracy of concentration profiles, we have compared our numerical code with the concentration and current profiles of a Glucose sensor [13]. Figure 4a shows glucose and H 2 O 2 concentration profiles in the multi-layer sensor system consisting of GO x /PPD and HAs/Fe 3+ as the first and second layers. Additionally, Figure 4b, depicts the dependency of the current density of the biosensor on the concentration of the glucose. All of the simulation parameters for comparison can be found in Table 2 of [13]. In addition, to ensure the outcomes in the presence of a perforated membrane and also the perforation level (α) and the level of the enzyme filling (γ) for a cylindrical geometry (Figure 2a), the results of the present study have been validated with the work of Baronas [20] in Figure 5. Finally, to compare the presented model with an experimental study, we have developed our FEM code to take into account the influence of the mediator layer and compared our modeling with the experimental data provided by Šimelevičius et al. [31] for mediator oxidization in an amperometric glucose biosensor. The maximum value of the residual sum of squares (RSS) for

Results and Discussion
To study the impact of geometry on the biosensor performance, seven geometries-including cylindrical, upward circular cone, downward circular cone, upward paraboloid, downward paraboloid, upward concave paraboloid, and downward concave paraboloid-are considered. The effects of the enzyme filling level and the perforation level in these holes on the transient and steadystate current and the half-time response of the biosensor have been investigated for a range of geometrical parameters. The values of parameters are considered as  Figure 7a,b for various geometries. It is clear that the geometry has a significant influence on the output current of the biosensor, and the transient current is higher for the upward geometries than the downward ones. In addition to this, increasing the perforation level boosts the output current strongly. However, the

Results and Discussion
To study the impact of geometry on the biosensor performance, seven geometries-including cylindrical, upward circular cone, downward circular cone, upward paraboloid, downward paraboloid, upward concave paraboloid, and downward concave paraboloid-are considered.  Figure 7a,b for various geometries. It is clear that the geometry has a significant influence on the output current of the biosensor, and the transient current is higher for the upward geometries than the downward ones. In addition to this, increasing the perforation level boosts the output current strongly. However, the

Results and Discussion
To study the impact of geometry on the biosensor performance, seven geometries-including cylindrical, upward circular cone, downward circular cone, upward paraboloid, downward paraboloid, upward concave paraboloid, and downward concave paraboloid-are considered. The effects of the enzyme filling level and the perforation level in these holes on the transient and steady-state current and the half-time response of the biosensor have been investigated for a range of geometrical parameters. The values of parameters are considered as a 1 = 10a 3 = S 0 = 1, D 3 = 2D 2 = 6, b 5 = 4b 2 = 16, b 4 = 7b 1 = 14 and σ 2 = 3.33 × 10 4 . The effects of the perforation level on the transient current are shown in Figure 7a,b for various geometries. It is clear that the geometry has a significant influence on the output current of the biosensor, and the transient current is higher for the upward geometries than the downward ones. In addition to this, increasing the perforation level boosts the output current strongly. However, the difference between the current in the upward and downward geometries decreased by increasing the level of perforation for the concave paraboloid and circular cone geometries. The maximum current can be obtained from the upward concave paraboloid holes. Figure 7c,d compares the effect of the enzyme filling level on the transient current of the biosensors with different geometries. It is clear that the γ affects both the biosensor response time and current. The current of the biosensor for high values of γ is up to three times higher than that obtained when the enzyme filling level approaches zero. Moreover, the difference between the upward and downward geometries vanishes for a low level of enzyme. Finally, the minimum current is produced in the cylindrical geometry.
difference between the current in the upward and downward geometries decreased by increasing the level of perforation for the concave paraboloid and circular cone geometries. The maximum current can be obtained from the upward concave paraboloid holes. Figure 7c,d compares the effect of the enzyme filling level on the transient current of the biosensors with different geometries. It is clear that the γ affects both the biosensor response time and current. The current of the biosensor for high values of γ is up to three times higher than that obtained when the enzyme filling level approaches zero. Moreover, the difference between the upward and downward geometries vanishes for a low level of enzyme. Finally, the minimum current is produced in the cylindrical geometry. Figure 8 depicts the effect of the perforation level and enzyme filling level on the nondimensional output current of the biosensor. Obviously, the current increases with the level of perforatio  Figure 8 depicts the effect of the perforation level and enzyme filling level on the non-dimensional output current of the biosensor. Obviously, the current increases with the level of perforation exponentially, and the current is highly sensitive to α when it approaches 1. Moreover, increasing the perforation level reduces the difference between the maximum and minimum steady-state currents of the hole-geometries, indicating that the impact of the hole-geometry on the output current decreases as α increases. In contrast with α, the output current increases with the increment of the enzyme filling level almost linearly.  Table 1 compares the effect of the perforation level and enzyme filling level on the half-time response of the biosensor. It is evident that the cylindrical and upward concave paraboloid holes require the lowest and highest time intervals to reach the steady-state condition. In addition to this, T0.5 is a monotonous function of the enzyme filling level and perforation level, and it increases with the increment of α and γ.   Table 1 compares the effect of the perforation level and enzyme filling level on the half-time response of the biosensor. It is evident that the cylindrical and upward concave paraboloid holes require the lowest and highest time intervals to reach the steady-state condition. In addition to this, T 0.5 is a monotonous function of the enzyme filling level and perforation level, and it increases with the increment of α and γ.

Conclusions
The influence of the perforated membrane geometry on the performance of a biosensor is studied in the case of a 2-D axisymmetric model of an amperometric biosensor. For this aim, seven practical geometries-including cylindrical, upward circular cone, downward circular cone, upward paraboloid, downward paraboloid, upward concave paraboloid, and downward concave paraboloid-are taken into consideration. In addition to this, the effects of the perforation level of the perforated membrane and the filling level of the enzyme level on the transient and steady-state current of the biosensor and the half-time response are studied. The results show that the geometry of the biosensor strongly influences the sensitivity and time response of the biosensor, specifying the potential of future experimental studies. The outcomes reveal that the maximum and minimum current can be obtained from the upward concave paraboloid and circular holes, respectively. Moreover, increasing the perforation level and enzyme filling level boosts the output current strongly. On the other hand, the time response of the studied biosensor for the cylindrical geometry is lower than that for other geometries.